System and method for time correlated multi-photon counting measurements

ABSTRACT

The invention provides a method and a measurement system for characterisation of luminescence properties, the method comprises irradiating the luminescent material with a pulse of excitation light, providing a triggering signal correlated to the pulse of excitation light; detecting with a photodetector such as a photomultiplier tube (PMT) a plurality of photons emitted from the luminescent material as result of the pulse of excitation light, the photodetector providing an output signal upon the event of detection of a photon; determining for each detected photon a photon arrival time and providing an output suitable for inputting to an analysing module wherein an output comprises zero, one, or more photon arrival time for each excitation, receiving said outputs in an analysing module; and determining. in the analysing module, characteristics properties of the luminescent material by performing a statistical analysis based on Bayesian inference.

The present invention relates to a method and system for measuringluminescence lifetime. More particularly, the invention relates to amethod and apparatus that significantly reduces the measurement time ina Time Correlated Multi-Photon Counting (TCMPC)-measurement.

BACKGROUND OF THE INVENTION

Measurement of the luminescence emission of substances has becomeincreasingly important as a sensitive analytical technique in a widerange of disciplines such as chemistry, biology, engineering, materialsscience and application areas ranging from pollution monitoring andproduction control in a chemical plant to the sequencing of DNA. In manyof these areas, distinguishing the weak but characteristic emission ofthe substance of interest from background light is difficult. Measuringthe lifetime of the light emission of the substance of interest canoften be used to discriminate the desired signal from a background. Inthis context we define the luminescence lifetime of a substance to bethe time taken for the intensity of light emitted in response to anexciting pulse of light to fall to 1/e (˜1/2.7) of its initial value.This characteristic time can range from less than one nanosecond tominutes, depending on the processes involved. Such processes can bedescribed as fluorescence, luminescence or phosphorescence, but if thedesired light emission has a lifetime which is different from light ofanother source, its detection can be enhanced.

One example of fluorescence measurements is in DNA sequencing, wheremolecular fragments are ‘labelled’ with a fluorescent dye and thefluorescence emission of a dye-labelled fragment is used to identifywhich amino acid base begins or ends a given sequence of amino acids.Another application is in the screening of potential drugs against atarget molecule. Typically changes in fluorescence intensity are used toidentify ‘hits’ in response to a chemical test or ‘assay’ of theinteraction between a potential drug and a target molecule. In thiscase, either the drug or the target or both may be labelled. In bothsequencing and screening, the characteristic lifetime of a fluorescentmolecule can be used to improve its detection.

A fluorescent material is a material which can be excited to a higherelectronic state and on relaxing towards its ground state emits aphoton. The time elapsed from the excitation until the emission of aphoton, the intensity of emitted photons, and the energy (frequency) ofthe photons are all properties that can be measured and used tocharacterize the material. The excited state is typically achieved bysubjecting the material to a short light pulse. The emission of photonsis a Poisson process and the intensity of the emitted photons iscontinuously decreasing, and typically follows an exponential decay. Auseful way of characterising a fluorescent material is to measure thetime elapsed between an excitation and the arrival of photons at adetector. λ(t) denotes the Poisson arrival rate at time t at all timesand is given by:λ(t)=Ae ^(−at)  equation 1where A denotes the intensity of a single fluorescent material, measuredin photons per second, at time 0 (the time of the excitation), alsoreferred to as the initial intensity, α denotes the decay constant ofthe fluorescence, measured in nepers/second, and 1/α is the fluorescencelifetime, which is often used to characterise a fluorescent material. Ifa sample has more than one fluorescent material present, the arrivalrate of photons can typically be described by a multiexponentialexpression.

A widely used technique of determining the fluorescence lifetime of amaterial is Time Correlated Single Photon Counting (TCSPC). In TCSPC asample is typically subjected to a short light pulse, which excites thematerial in the sample. The material will, if fluorescent, emit a photona short time after the excitation. The emitted photon is detected andthe time that has elapsed between the excitation light pulse and thearrival of the first emitted fluorescence photon is recorded by themeasurement system. The technique is only capable of counting a maximumof one photon per excitation. In practice less than one photon perexcitation are counted since, to get a representative distribution ofarrival times, a low intensity of the exciting light pulse is used,giving a low probability of photon emission. Therefore, to get areliable measure of the fluorescence lifetime characteristic of aparticular material, using TCSPC, i.e. to acquire a fluorescence decaycurve, the process of exciting the material in the sample and detectingthe first emitted photon, has to be repeated a large number of times. Ameasurement can typically involve several tens of thousandsexcitation-photon detection cycles. The TCSPC method has beencontinuously improved, for example U.S. Pat. No. 5,990,484 teaches amethod and apparatus for counting a total number of emitted photons,allowing a higher intensity of the exciting light pulse, and U.S. Pat.No. 4,686,371 teaches a method of compensating for variations in thelength and intensity of the exciting light pulse.

As above indicated, TCSPC can, even if the mentioned improvements areutilised, be a time consuming measurement. It is realised in the artthat the method is wasteful in that only the first emitted photon isdetected and used for the determination of the fluorescence lifetime,even though the subsequently emitted photons carry information usefulfor the determination of the characteristic lifetime. A number ofapproaches have been suggested to extend the TCSPC to also detect andrecord subsequent photons. The most recent development of theTCSPC-technique includes various multiplexing techniques, in which, inprinciple, a plurality of detectors are connected to one or moreanalysing means, examples described in Multiplexed single-photoncounting by Scubling et al., Rev. Sci. Instrum. 67 (6), 2238-2246 June1996. Alternatively multi-element detectors have been suggested, forexample by Browne et al. in A 100 ns anti-coincidence circuit formulti-channel counting systems, Meas. Sci. Technol. 6 (6), 1487-1491,June 1996. The suggested Multiplexed Time Correlated Single-PhotonCounting (M-TCSPC)—techniques represent significant improvements overthe TCSPC—techniques, but have the drawback that the number of photonsthat it is possible to detect following one excitation pulse is limitedby the number of detectors (or detector elements) used; therefore alarge number of excitation cycles is still needed to determine thefluorescence lifetime of the material. In addition multielementdetectors are costly, as are systems providing a large number of singledetectors in a multiplexing arrangement

The known apparatus and measurement techniques thus suffer, as abovedescribed, from a number of disadvantages and limitations, the mostimportant being the long measurement time. Thus, there is a need in theart to provide a measurement system and method that shortens measurementtime.

SUMMARY OF THE INVENTION

The object of the present invention is to provide a system and a methodthat significantly reduces the measurement time in a major part of theapplications. The inventive system is defined in claim 10 and the methodand corresponding software module in claims 1 and 9, respectively.

One advantage afforded by the present invention is to be able to use alimited dataset by using a Bayesian algorithm to make better use of theavailable data.

Yet another advantage of the invention is to be able to take intoaccount the PMT dead time in the lifetime analysis or any otherinstrumental parameters, such as e.g non-delta illumination, etc.

Yet another advantage of the invention is the possibility, as a resultof a measurement, of giving the posterior probability distribution of afluorescence lifetime and/or intensity in contrast to the prior art,which typically only gives a single value as a measure of the lifetimeand/or intensity.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention will now be described in detail with reference to thedrawing figures, in which

FIG. 1 is a schematic illustration of a single channel TCMPC systemaccording to the invention;

FIG. 2 is a schematic illustration of a multiple channel TCMPC systemaccording to a preferred embodiment of the invention;

FIG. 3 shows a flowcharts of the measurement and analysing methodaccording to the invention;

FIG. 4 a is a graph showing an example of a trace from the photon timingunit;

FIG. 4 b is a graph showing an example of a trace from the photon timingunit from the multi-channel TCMPC system;

FIG. 5 is a graph showing an example of three prior distributionsentered by a user;

FIG. 6 is a graph showing an exemplary output from algorithms A3-A8.

DETAILED DESCRIPTION OF THE INVENTION

The inventive measurement system provided by the present invention willbe described with reference to the schematic illustrations of FIG. 1-2.FIG. 1 shows a single channel TCMPC system according to the invention.Fluorescent material is contained in a sample 110. An excitation lightsource 120 for repeatedly irradiating the sample with pulses ofexcitation light is provided. The light source preferably a pulsedlaser. A light source driving and controlling unit 130 comprises thelaser power supply and is arranged to generate a triggering signal 140to a photon timing unit 150. The repetition rate of the excitation lightshould be sufficiently low to allow a total decay of the fluorescence ofthe sample before the next exciting pulse. A pulsing frequency of forexample around 200 kHz would commonly be appropriate. The duration ofthe pulse of excitation light should be significantly shorter than thefluorescence lifetime of the fluorophore in order to obtain reliablelifetime measurements. A typical fluorophore has a fluorescence lifetimein the range of 1-100 ns and a suitable length of the light pulse wouldbe in the order of 0.5 ns. The frequency of the pulse is an example of aparameters typically adjusted for the sample under investigation (thelength of the light pulse is an intrinsic characteristic of the laserand typically cannot be easily changed). In addition, the light source120 and the light source driving and controlling unit 130, shouldpreferably have the facility to adjust other parameters such as thenumber of pulses and the intensity of the excitation light in order toaccount for e.g. the amount of fluorophore in the sample, the neededaccuracy in the result etc. Light sources, e.g. pulsed lasers, with thecharacteristics described above and with driving and controlling unitstherefore are known in the art and commercially available; for example,the NanoLED-07 picosecond laser source from IBH.

Positioned adjacent to a sample chamber containing the sample 110 is aphotodetector 160, which has the purpose of detecting the emittedfluorescent photons. Different optical components are placed between thesample and the detector. For example several lenses are used to maximisethe amount of fluorescent light collected from the sample and to focusthe light onto the detector. Furthermore dichroic mirrors and filtersare used to select a range of wavelength and prevent the excitationlight from reaching the detector and only let the light having a rangeof wavelength corresponding to the fluorescence spectrum of thefluorophore reach the detector. An optical component is also used tosplit the beam of the fluorescent light into several beams that are thendirected to the different detectors. This splitting can be achieved indifferent ways, e.g. by using a beam splitter cube or by usingmulti-branch fibre optics. These are well known in the art and arecommercially available (for example, from Oriel Instrument, 150 LongBeach Boulevard, Stratford, Conn., USA or from Melles Griot, St ThomasPlace, Ely, Cambridgeshire, England). The photodetector preferablycomprises a photon counting photomultiplier tube (PMT) but otherdetectors such as an avalanche photodiode can be used. The signal fromthe photodetector is typically amplified in a pre-amplifier 170. Afteramplification, the signal can go through a discriminator unit thatremove the unwanted noise from the signal and only leave the electricalpulses generated by the photons on the PMT. The discriminator unit isthen connected to the photon timing unit 150. The photon timing unitcomprises for example a very fast analogue-to-digital converter (A/Dconverter) and a memory for storing datapoints. As discussed above, itshould be possible for more than one photon per excitation, resultingfrom the excitation light pulse, to be recorded by the photodetector 160and the photon timing unit 150, not only the first emitted photon as inthe TCSPC based techniques. This put certain requirements on thephotodetector 160 and the photon timing unit 150; A photon detected bythe PMT will give rise to an output pulse that has a duration of c.a.4-6 ns with the current technology. To identify a PMT pulse position intime at least 4-8 data points per pulse are needed, in accordance withthe method described further below. In order to detect a plurality ofsuch pulses the A/D converter needs to have a sampling frequency atleast of around 2 GS/s and the ability to store data pointscorresponding to an acquisition period of approximately 100 ns. A/Dconverters with such extreme capabilities have only recently been madecommercially available. An example of a suitable A/D converter is theAcquiris DP210 (2 GS/s) from Acquiris Europe, 18 chemin des Aulx, 1228Plan-les-Ouates, Geneva, Switzerland. The collected data points areanalysed by an arrival time determination module 155, preferablyrealised as a software program module residing either within the photontiming unit 150, or alternatively within an external computer 180, todetermine the arrival time of the fluorescence photons. Examples ofsuitable algorithms for arrival time determination will be given below.Arranged to receive and analyse a dataset of photon arrival times fromthe arrival time determination module, is an analysing module 185, whichpreferably is a software program package implemented in and executed onthe computer 180. The computer can be a general purpose computer such asa PC or a workstation or a special purpose machine. The computer 180 ispreferably also used for controlling the measurements, monitoring theequipment etc. and therefore is equipped with communication means (notshown) connected to other units of the measurement system. Installed inand executing on the computer, is then also a software program modulefor measurement control 190. As appreciated by the skilled in the art,the installation of and execution of the software modules 155, 185 and190 can be done in various ways and in various computer configurations.For example the measurement control unit 190 can be executed in a PCsuitable for laboratory conditions, the arrival time determinationmodule 155 can be incorporated within the photon timing unit, and theanalysing module 185 in a special purpose machine designed for the highcomputational speed needed for certain analysis. The inventive methodsutilised by the analysing module 185 will be described in detail below.

Photodetectors such as PMTs give a signal that typically last 4-5 ns.During that time the PMT is incapable of detecting any new photons. Thisis often referred to as PMT dead time. If photons arrive during the PMTdead time, i.e. less than 5 ns apart, the measurement system will not beable to record them. As will be described further below, the analysingmethod can be made to take into account the PMT dead time, but the factthat the system is incapable of recording photons during certain timeintervals will prolong the needed measurement time.

In one embodiment of the present invention the measurement system hasbeen provided with a plurality of photodetectors (PMT) and amplifiersconnected to them, giving a multi-channel TCMPC. In FIG. 2 this isexemplified with a dual-channel system, having a first PMT 200, and asecond PMT 205, connected to preamplifiers 210 and 215, respectively.The preamplifiers are connected to a common photon timing unit 150 withindependent channels. Alternatively a plurality of photon timing unitscan be used. If one of the PMT is recovering after detecting a photon,the other may still be capable of detecting a subsequent photon. Theeffective dead time of the system is thus reduced. More than twoPMT—amplifier channels can be used to further reduce the effective PMTdead time.

An experimental procedure using the inventive methods and measurementsystem according to the invention will be described with references tothe flowchart of FIG. 3. Variations of the experimental procedure arepossible and the procedure now described should be considered as one ofmany embodiments possible. The aim of the measurement is typically todetermine the initial intensity of fluorescence, the fluorescencelifetime (1/α) or the posterior probability distribution of theintensity or lifetime of one or a plurality of fluorophores contained inthe sample 110.

The experimental procedure comprises:

-   -   300: In a first step 300, the user enters, to the measurement        controlling module 190 and the analysing module 185, a number of        parameters used for controlling the measurement, choosing an        appropriate analysis algorithm depending primarily on prior        knowledge of the sample, and parameters describing the        performance of the equipment such as PMT dead time. Parameters        important for the inventive methods described below comprise:        -   N—the number of excitations (laser pulses);        -   α_(i)—the reciprocal of the fluorescence lifetimes (if            known) and/or the prior distribution of lifetimes, P(α_(i));        -   A_(i)—the fluorescence initial intensities (if known) and/or            the prior distribution of intensities, P(A_(i));        -   δ—the PMT dead time;        -   τ—the time the detector is blocked after an excitation;        -   υ—the duration of recording after each excitation        -   A1, A2, A3, A4, A5, A6, A7, A8—Analysis algorithms        -   Many of the parameters are not changed for each measurement            and default values are preferably stored in a database.    -   305: The measurement starts by illuminating the sample 110 with        a pulsed laser 120.    -   310: In a next step, 310, the emitted photons are collected by        PMT tube(s) 160 and the signal is amplified by a pre-amplifier        170 and possibly the signal is discriminated to remove unwanted        noise.    -   315: In step 315, the arrival time of the emitted photons are        determined by the photon timing unit 150 and the arrival time        determination module 155, with one of the methods that will be        described below. The arrival time determination module 155        outputs for each excitation n a dataset D_(n), comprising the        photon arrival times, to the analyzing module 185.    -   320: In step 320, if the predetermined number of excitations N        is reached the procedure goes to step 325, if not, steps 300-320        are repeated.    -   325: In step 325, the datasels D_(n) are analysed with one (or        more) of the algorithms A1-A8, the algorithms being based on        Bayesian inference.    -   330: In a final step 330, the result is presented to the user.        Depending on the algorithm chosen and the preferences of the        user the output can be distributions of the fluorophore        lifetimes or initial intensities, joint distributions, a measure        of a fluorophore lifetime etc.

Alternatively the measurement procedure can stop then a predefinedaccuracy is achieved. In this case the analysis using Bayesian inference(step 325) is performed within the main loop of the algorithm. If, forexample, the measure of a fluorescence lifetime has the predefinedaccuracy, the loop is halted. The procedure can of course in additionhave a limit to the number of pulses used, which terminates themeasurement even if the specified accuracy has not been met.

Arrival Time Determination

The method of determining photon arrival times performed in the arrivaltime determination module 155, will be described with reference to FIG.4 a-b. FIG. 4 a gives an example of an output, commonly known as atrace, from the photon timing unit 150. The initial broad peak in thetop curve corresponds to the triggering signal 140 sent from the lightsource driving and controlling unit 130 to the photon timing unit 150,when the excitation light source 110 irradiates the sample. It should benoted that this triggering signal itself is much longer (in this case itlasts for about 50 ns) than the actual laser pulse (that last for only0.5 ns). Only the rising part of the signal is used by the photon timingunit 150 to determine the trigger to start the time base. The rest ofthe triggering signal is ignored by the photon timing unit. It must benoted that the length of the triggering signal does not correspond tothe length of the laser pulse. If the laser pulse is longer than theresolution of the photon timing unit, one can block the photodetectorduring a time τ. An emitted fluorescence photon is detected by the PMTand give rise to the pulse from the PMT, amplified in the pre-amplifier(solid line on FIG. 3), discriminated to remove the noise (dotted lineon FIG. 3) and digitised in the A/D-converter, to show as the peak 210.It should be noted that the discriminator introduces a delay in thesignal, explaining the shift between the position of the pulse in thesolid and dotted line. This is not a problem, as this time delay isconstant, it can be subtracted subsequently by the timing unit The timeelapsed between the triggering signal and the appearance of the firstpeak 210 is the arrival time of the first fluorescence photon. Thesecond peak 220 corresponds to a subsequent fluorescence photon.

A method of determining the arrival time of the fluorescence photonswill now be described with references to FIG. 4 a-b The trace is a listof data points, each point representing an output voltage from the PMTat a given time. In the embodiment with a plurality of photodetectors acorresponding plurality of parallel traces will result. In a first step505 a zero time is defined as the beginning of the triggering pulse (forexample, at the time when the trigger pulse reaches 50% of its maximumamplitude). The second step 510 is discrimination by selecting a part ofthe trace that is larger than a threshold value (this can be done bysoftware or by hardware in the discriminator). In the following step 515a local maximum of the voltage is defined as a point that is larger thanboth the point directly preceding it and the point directly followingit. The time corresponding to the local maximum is recorded in step 520.This is defined as the arrival time of the photon. The method ispreferably realised by a software program executed either in the photontiming unit 150 or in an external computer 180. The result is a datasetD_(n) of measured arrival times for each excitation pulse n. Asappreciated by the skilled in the art other methods of finding arrivaltimes in a trace of datapoints can be utilised. For example, a firstpoint larger than a preset threshold, or the point with the largestslope (largest derivative) in the peak, can be used to define an arrivaltime. The method chosen is typically a trade-off of the requiredprecision in the definition of the arrival time and the computationaltime needed to determine the arrival time.

FIG. 4 b shows the two traces, and the trigger signal from the laserpulse, resulting from the measurement system of FIG. 2. Both themulti-channel and the single-channel TCMPC can further be provided withdiscriminators, one for each channel, or other noise reduction means.

It should be noted that with all these methods, there is sometime anartefact. The algorithm may identify two points (or more) in a singlepulse that satisfy the selection criteria. In order to remove theartefactual points, a subroutine would check whether several selectedpoints are spaced with a time less than the PMT dead time and onlyselect the first point in the group.

Bayesian Inference

The method of the present invention, which will be illustrated by theembodiments, algorithms A1-A8, are based on the principle of Bayesianinference. Bayesian inference offers a way of determining the(posterior) probability distribution of a variable (or vector ofvariables) in the light of some data taking into account prior knowledgeof the probability of the variable taking various values. As an example,let us suppose that θ represents what we want to know. Usually θ will bea vector of many numbers; it might for example be the vector ofintensities of a number of fluorophores. Let us suppose that Drepresents what we do know. D might for example the arrival times ofsome photons. Now, we know that for any specific value of θ, whichspecifies the vector of intensities, that some values of D are morelikely and others less likely, and because we know the manner in whichthe data arises from θ, we can say exactly how much more or less likely.In other words, we can evaluate, for any specific values of θ and D, thevalue of P(D|θ), the probability of D given θ. We also know, before weeven receive D, that some values of θ are more likely than others. Forthe photon counting example negative intensities can be excluded, as canvalues of the intensities that are extremely large (e.g. >1 e15 persecond) or too small. This can all be precisely specified as a priordistribution on θ, denoted P(θ).

What we actually want to know is the value of θ, given some specificvalue of D. This is given by Bayes' theorem: $\begin{matrix}{{P\left( \theta \middle| D \right)} = \frac{{P(\theta)}{P\left( D \middle| \theta \right)}}{\int{{P(\theta)}{P\left( D \middle| \theta \right)}{d\theta}}}} & {{equation}\quad 2}\end{matrix}$which tells us that given that particular value of D, how likely eachvalue of θ is. In eq. 2 P(θ) is referred to as the Prior, P(D|θ) as theLikelihood, P(θ|D) as the Posterior, and ∫P(θ)P(D|θ)dθ as the Evidence(for the model).

A posterior distribution derived by Bayesian inference may by used indifferent ways: View the entire posterior distribution—possible in somecases (dimensionality not more than 2) and impossible in others; MAP(maximum a posteriori)—find the value of θ where the posteriorprobability density is greatest (if possible); and report just thatvalue of θ, Mean posterior value—report the mean posterior value of θ,θθP(θ|D)dθ; Random samples—drawing random samples from the distributionP(θ|D). The method of random samples is often the preferred way of usingthe Posterior. Random samples can be independent or non-independent.

The Bayesian principles are utilized in the inventive method of thepresent invention for substantially shortening the measurement timeneeded for e.g. determination of a fluorescence lifetime. The reductionin measurement times arises from the fact that by using Bayesianinference adapted to the application of TCMPC fewer photons are neededto characterize a fluorophore compared to the prior art methods. Themethod fully utilizes the advantage afforded by the inventive system ofdetecting a plurality of photons for each excitation pulse. Depending onthe prerequisite of the measurement, e.g. the number of differentfluorophores present in the sample, different embodiments of the methodwill be appropriate, corresponding to the mentioned algorithms A1-A8.The algorithms are preferably comprised in the analysing module 185, assoftware programs implemented in and executed on the computer 180(although it would be possible to embody them in other ways). In themeasurement procedure outlined with references to FIG. 3, the programsbased on the algorithms are addressed in the initialisation step 300 andthe analysis step 325. A brief theoretical background will be given toeach algorithm and the algorithms will be presented as pseudo-code withreferences to accompanying flowcharts.

Algorithm A1, Unknown Fluorophore Lifetimes, Single Fluorophore, NoDeadtimes.

The decay rate of a single fluorophore is given by equation 1. To useBayesian inference, we must have a prior distribution on the unknowns.We will use Gamma distributions on both A and α, as this is bothsufficiently expressive and results in relatively simple mathematics.The Gamma distribution has two parameters, m, the shape parameter, andr, the scale parameter, both of which must be positive. Its probabilitydensity is given by $\begin{matrix}{{\Gamma_{m,r}(x)} = {{\frac{r^{m}x^{m - 1}e^{- {rx}}}{\Gamma(m)}{or}\quad{for}\quad\log\quad{x:{\Gamma_{m,r}\left( {\log(x)} \right)}}} = \frac{r^{m}x^{m}e^{- {rx}}}{\Gamma(m)}}} & {{Equation}\quad 3}\end{matrix}$

Varying r scales the whole distribution; varying m changes the shape ofthe distribution, allowing very broad distributions when m is small andvery narrow ones when m is big. We can therefore choose appropriateparameters, and set $\begin{matrix}{{{P(A)} = \frac{r_{A}^{m_{A}}A^{m_{A} - 1}e^{{- r_{A}}A}}{\Gamma\left( m_{A} \right)}}{and}} & {{Equation}\quad 4} \\{{P(a)} = {\frac{r_{a}^{m_{a}}a^{m_{a} - 1}e^{{- r_{a}}a}}{\Gamma\left( m_{a} \right)}.}} & {{Equation}\quad 5}\end{matrix}$

As well as a prior, we must be able to calculate the likelihood of theobserved data. Let us suppose that during a excitation n, photons arriveat times t_(n,1),t_(n,2), . . . ,t_(n,K) _(n) , and at no other times.Then we have the likelihood on the data D_(n) from excitation n given by$\begin{matrix}{{{P\left( {\left. D_{n} \middle| A \right.,a} \right)} = {\left( {\prod\limits_{k = 1}^{K_{n}}\quad\left( {Ae}^{- {at}_{n,k}} \right)} \right){P\left( {\left. {{no}\quad{other}\quad{photons}} \middle| A \right.,a} \right)}}},} & {{Equation}\quad 6}\end{matrix}$since in a Poisson process all events occur independently of each other.In a Poisson process of constant rate λ, the probability of no photonsoccurring in the interval from t₁, to t₂ is e^(−λ(t) ² ^(−t) ¹ ⁾. Wherethe rate is not constant, this probability is e $\begin{matrix}\begin{matrix}{{P\left( {\left. D_{n} \middle| A \right.,a} \right)} = {\left( {\prod\limits_{k = 1}^{K_{n}}\quad\left( {Ae}^{- {at}_{n,k}} \right)} \right)e^{- {\int_{0}^{\infty}{{\lambda{(t)}}{dt}}}}}} \\{= {A^{K_{n}}e^{{{- a}{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}} - {\int_{0}^{\infty}{{Ae}^{- {at}}{dt}}}}}} \\{= {A^{K_{n}}e^{{{- a}{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}} - \frac{A}{a}}}}\end{matrix} & {{Equation}\quad 7}\end{matrix}$

We therefore have e^(−∫_(t₁)^(t₂)λ(t)dt).and combining all the data from N excitations together (independent fromexcitation to excitation), we have${P\left( {\left. D \middle| A \right.,a} \right)} = {{A^{K}e^{{- {a{({\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}})}}} - \frac{NA}{a}}{where}\quad K} = {\sum\limits_{n = 1}^{N}\quad{K_{n}.}}}$

Now, applying Bayes' theorem to determine the posterior, we find that$\begin{matrix}\begin{matrix}{{P\left( {A,\left. a \middle| D \right.} \right)} = \frac{{P\left( {A,a} \right)}{P\left( {\left. D \middle| A \right.,a} \right)}}{\int{{P\left( {A,a} \right)}{P\left( {\left. D \middle| A \right.,a} \right)}{dAda}}}} \\{\propto {\frac{r_{A}^{m_{A}}A^{m_{A} - 1}e^{{- r_{A}}A}}{\Gamma\left( m_{A} \right)}\frac{r_{a}^{m_{a}}a^{m_{a} - 1}e^{{- r_{a}}a}}{\Gamma\left( m_{a} \right)}A^{K}e^{{- {a{({\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}})}}} - \frac{NA}{a}}}} \\{\propto {A^{m_{A} - 1}e^{{- r_{A}}A}a^{m_{a} - 1}e^{{- r_{a}}a}A^{K}e^{{- {a{({\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}})}}} - \frac{NA}{a}}}} \\{= {A^{K + m_{A} - 1}a^{m_{a} - 1}e^{- {A{({r_{A} + \frac{N}{a}})}}}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}}}\end{matrix} & {{Equation}\quad 8}\end{matrix}$which gives the joint posterior distribution and hence $\begin{matrix}\begin{matrix}{{P\left( a \middle| D \right)} = {\int{{P\left( {A,\left. a \middle| D \right.} \right)}{dA}}}} \\{\propto {\int{A^{K + m_{A} - 1}a^{m_{a} - 1}e^{- {A{({r_{A} + \frac{N}{a}})}}}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}{dA}}}} \\{= {a^{m_{a} - 1}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}{\int{A^{K + m_{A} - 1}e^{- {A{({r_{A} + \frac{N}{a}})}}}{dA}}}}} \\{= \frac{a^{m_{a} - 1}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}{\Gamma\left( {K + m_{A}} \right)}}{\left( {r_{A} + \frac{N}{a}} \right)^{K + m_{A}}}}\end{matrix} & {{Equation}\quad 9}\end{matrix}$giving the marginal distribution of a (the reciprocal of the fluorophorelifetime), while $\begin{matrix}{{P\left( {\left. A \middle| a \right.,D} \right)} = {\frac{\left( {r_{A} + \frac{N}{a}} \right)^{K + m_{A}}}{\Gamma\left( {K + m_{A}} \right)}A^{K + m_{A} - 1}{e^{- {A{({r_{A} + \frac{N}{a}})}}}.}}} & {{Equation}\quad 10}\end{matrix}$gives the conditional distribution of initial intensity A.

With the theoretical background above, and with reference to theflowchart of FIG. 3, an algorithm suitable for software implementationwill be described. The algorithm will be denoted A1. As previousdiscussed, although the lifetime 1/α and the initial intensity A of thefluorophore is unknown, the user normally has some vague knowledge ofthe possible range of these values, that can be used to generate theprior distribution needed for the Bayesian inference.

Algorithm A1:

-   500: In a first initialising step 500, which is included in step 300    of the flowchart of FIG. 3, the user enters functions, values or    distributions describing α and A. If possible, the program may plot    the result, and the user acknowledge the outcome.-   505: In a measuring step 505, which corresponds to the steps    305-320, datasets of arrival times are collected D_(n).-   510: The Bayesian inference analysis, step 510 (corresponding    to 325) comprises the substeps (510: 1-5) of:-   510: 1 Set grid of A, α (or log A, log α).-   510: 2 Calculate $\begin{matrix}    {{P\left( {A,\left. a \middle| D \right.} \right)} = {A^{K + m_{A} - 1}a^{m_{a} - 1}e^{- {A{({r_{A} + \frac{N}{a}})}}}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}}} & \left( {{eq}.\quad 8} \right)    \end{matrix}$    for each combination of values of A and α-   510: 3 Normalize the joint distribution P(A,α|D) by:    -   multiplying each grid point value of P(A,α|D) by the product        measure appropriate to that point,    -   summing over all axes, dividing the original distribution by the        result of the summation.-   510: 4 Marginalize P(A,α|D) to get P(α|D) and P(A|D) (eq. 9,10) by,    for each margin:    -   multiplying each grid point value of P(A,α|D) by the product        measure appropriate to that point,    -   summing over all axes except the axis whose margin is wanted,    -   dividing each grid point value of the resulting distribution by        the marginal measure appropriate to that point-   520: In an outputting step 520 (corresponding to 330) the resulting    distributions are presented to the user as graphs, for example.    Other measures can be calculated from the distribution for example    the marginal means with the 5^(th) and 95^(th) centiles, an often    quoted measure of the lifetime of a fluorophore.

The computational steps above usually actually need to be done byworking in logp or −logp rather than in P, in order to avoid underflowand overflow. The skilled in the art will appreciate that in thecalculating steps a number of obvious computational savings may be made,such as avoiding recalculating parts not changing during the summations.

The normalising step 510: 3 may preferably be performed in the followingmanner:

The grid of values of A are x_(h)=(x₁,x₂, . . . ,x_(H)) and those of αare y_(j)=(y₁,y₂, . . . ,y_(j)) Associate a measure with each x_(h) andy_(j):${m\left( x_{h} \right)} = {{\frac{x_{h + 1} + x_{h}}{2} - {\frac{x_{h} + x_{h - 1}}{2}{and}\quad{m\left( y_{j} \right)}}} = {{\frac{y_{j + 1} + y_{j}}{2} - {\frac{y_{j} + y_{j - 1}}{2}{and}\quad{m\left( {x_{h},y_{j}} \right)}}} = {{m\left( x_{h} \right)}{m\left( y_{j} \right)}}}}$

The normalised joint distribution will be given by${P\left( {{A = x_{h}},{a = \left. y_{j} \middle| D \right.}} \right)}:=\frac{P\left( {{A = x_{h}},{a = \left. y_{j} \middle| D \right.}} \right)}{\sum\limits_{h}\quad{\sum\limits_{j}\quad{{P\left( {{A = x_{h}},{a = \left. y_{j} \middle| D \right.}} \right)}{m\left( {x_{h},y_{j}} \right)}}}}$where all the symbols indicate the values stored at the grid points andthe assignment symbol :=indicates new values to be stored in place ofthe old (which were only correct up to a constant of proportionality).

The marginalising step 510: 4 may preferably be performed bycalculating:${P\left( {A = \left. x_{h} \middle| D \right.} \right)}:={\frac{\sum\limits_{j}\quad{{P\left( {{A = x_{h}},{a = \left. y_{j} \middle| D \right.}} \right)}{m\left( {x_{h},y_{j}} \right)}}}{m\left( x_{h} \right)}{and}}$${P\left( {a = \left. y_{j} \middle| D \right.} \right)}:=\frac{\sum\limits_{j}\quad{{P\left( {{A = x_{h}},{a = \left. y_{j} \middle| D \right.}} \right)}{m\left( {x_{h},y_{j}} \right)}}}{m\left( y_{j} \right)}$

The optional calculation of marginal means with the 5^(th) and 95thcentiles in the outputting step 520 may preferably be performed bycalculating the cumulative density functions;${P\left( {A \leqq x_{h}} \right)} = {{{\sum\limits_{h^{\prime} = 1}^{h}\quad{{P\left( {A = x_{h^{\prime}}} \right)}{m\left( x_{h^{\prime}} \right)}\quad{and}\quad{find}\quad{greatest}\quad h\quad{such}\quad{that}\quad{P\left( {A \leqq x_{h}} \right)}}} \leqq {0.05\quad{and}\quad{P\left( {A \geqq x_{h}} \right)}}} = {{\sum\limits_{h^{\prime} = 1}^{h}\quad{{P\left( {A = x_{h^{\prime}}} \right)}{m\left( x_{h^{\prime}} \right)}\quad{and}\quad{find}\quad{least}\quad h\quad{such}\quad{that}\quad{P\left( {A \geqq x_{h}} \right)}}} \leqq 0.05}}$Single Fluorophore, with Photon-Induced, Initial, and Final Deadtimesand Unknown Lifetime

One of the major advantages with Bayesian inference is the possibilityof introducing parameters describing e.g. equipment performance andartifacts in the analysis. The result will be the true distributionstaking these parameters into account, rather than compensating a resultwith heuristic factors reflecting the artifacts, as in prior art.

Adopting the same notation as before, suppose that the detector is deadfor a duration δ after each photon is received, for a duration τ afterthe excitation, and at all times longer than υ after excitation. Again,we must calculate the likelihood of the observed data. Let us supposethat during a excitation n, photons arrive at times t_(n,1),t_(n,2), . .. ,t_(n,K) _(n) , and at no other times. As before we have thelikelihood on the data D_(n) from excitation n given by${{P\left( {\left. D_{n} \middle| A \right.,a} \right)} = {\left( {\prod\limits_{k = 1}^{K_{n}}\quad\left( {Ae}^{- {at}_{n,k}} \right)} \right){P\left( {\left. {{no}\quad{other}\quad{photons}} \middle| A \right.,a} \right)}}},$since in a Poisson process all events occur independently of each other.The difference between this expression and the corresponding one is inthe meaning of the phrase “no other photons”; in the preceding sectionit meant no other photons ever, while now it just means no other photonsin the time that the detector is live.

As before, this probability is e^(−∫_(T)λ(t)  𝕕t),where T now denotes the set of times at which the detector is live. Wetherefore have${P\left( {\left. {{no}\quad{other}\quad{photons}} \middle| A \right.,a} \right)} = {e^{- {({{\int_{r}^{t_{n,1}}{{\lambda{(t)}}\quad{\mathbb{d}t}}} + {\sum\limits_{k = 1}^{K_{n} - 1}\quad{\int_{t_{n,k}}^{t_{n,k} + 1}{{\lambda{(t)}}\quad{\mathbb{d}t}}}} + {\int_{t_{n,K_{n}}}^{v}{{\lambda{(t)}}\quad{\mathbb{d}t}}}})}} \approx e^{- {({{\int_{r}^{v}{{\lambda{(t)}}\quad{\mathbb{d}t}}} - {\sum\limits_{k = 1}^{K_{n}}\quad{\int_{t_{n,k}}^{t_{n,k} + \delta}{{\lambda{(t)}}\quad{\mathbb{d}t}}}}})}}}$where the last line follows if we make the assumption that none of thedeadtimes overlap with the period following υ. Thus $\begin{matrix}{{{P\left( {\left. {{no}\quad{other}\quad{photons}} \middle| A \right.,a} \right)} \approx e^{- {({{\int_{r}^{v}{{\lambda{(t)}}\quad{\mathbb{d}t}}} - {\sum\limits_{k = 1}^{K_{n}}\quad{\int_{t_{n,k}}^{t_{n,k} + \delta}{{\lambda{(t)}}\quad{\mathbb{d}t}}}}})}}} = e^{- {({{\int_{r}^{v}{{Ae}^{- {at}}\quad{\mathbb{d}t}}} - {\sum\limits_{k = 1}^{K_{n}}\quad{\int_{t_{n,k}}^{t_{n,k} + \delta}{{Ae}^{- {at}}\quad{\mathbb{d}t}}}}})}}} \\{= {e^{\frac{A}{a}{({{\lbrack e^{- {at}}\rbrack}_{r}^{v} - {\sum\limits_{k = 1}^{K_{n}}\quad{\lbrack e^{- {at}}\rbrack}_{n,k}^{n,k^{+ \delta}}}})}} = e^{{- \frac{A}{a}}{({{({e^{- {at}} - e^{- {av}}})} - {\sum\limits_{k = 1}^{K_{n}}\quad{({e^{- {at}_{n,k}} - e^{- {a{({t_{n,k} + \delta})}}}})}}})}}}} \\{= {e^{{- \frac{A}{a}}{({{({e^{- {ar}} - e^{- {av}}})} - {\sum\limits_{k = 1}^{K_{n}}\quad{(e^{- {{at}_{n,k}{({1 - e^{- {a\delta}}})}}})}}})}} = e^{{- \frac{A}{a}}{({{({e^{- {ar}} - e^{- {av}}})} - {{({1 - e^{- {a\delta}}})}{\sum\limits_{k = 1}^{K_{n}}\quad e^{- {at}_{n,k}}}}})}}}}\end{matrix}$so combining all the data from N excitations together (independent fromexcitation to excitation), we have $\begin{matrix}{{P\left( {\left. D \middle| A \right.,a} \right)} = {\prod\limits_{n = 1}^{N}\quad\left( {\left( {\prod\limits_{k = 1}^{K_{n}}\quad\left( {Ae}^{- {at}_{n,k}} \right)} \right)e^{{- \frac{A}{a}}{({{({e^{- {ar}} - e^{- {av}}})} - {{({1 - e^{- {a\delta}}})}{\sum\limits_{k = 1}^{K_{n}}\quad e^{- {at}_{n,k}}}}})}}} \right)}} \\{= {A^{K}e^{{- {a{({\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}})}}} - {\frac{A}{a}{({{N{({e^{- {ar}} - e^{- {av}}})}} - {{({1 - e^{- {a\delta}}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{- {at}_{n,k}}}}}})}}}}}\end{matrix}.$

Now, applying Bayes' theorem to determine the posterior, we find that$\begin{matrix}{{P\left( {A,{\alpha ❘D}} \right)} = {{\frac{{P\left( {A,\alpha} \right)}{P\left( {{D❘A},\alpha} \right)}}{\int{{P\left( {A,\alpha} \right)}{P\left( {{D❘A},\alpha} \right)}{\mathbb{d}A}\quad d\quad\alpha}} \propto {\frac{r_{A}^{m_{A}}A^{m_{A} - 1}e^{{- r_{A}}A}}{\Gamma\left( m_{A} \right)}\frac{r_{\alpha}^{m_{\alpha}}\alpha^{m_{\alpha} - 1}e^{{- r_{\alpha}}\alpha}}{\Gamma\left( m_{\alpha} \right)}A^{K}e^{{- {\alpha{({\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}t_{n,k}}})}}} - {\frac{A}{\alpha}{({{N{({{\mathbb{e}}^{- {\alpha\tau}} - {\mathbb{e}}^{- {\alpha\upsilon}}})}} - {{({1 - {\mathbb{e}}^{- {\alpha\delta}}})}{\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}{\mathbb{e}}^{{- \alpha}\quad t_{n,k}}}}}})}}}} \propto {A^{m_{A} - 1}e^{{- r_{A}}A}\alpha^{m_{\alpha} - 1}e^{{- r_{\alpha}}\alpha}A^{K}e^{{- {\alpha{({\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}t_{n,k}}})}}} - {\frac{A}{\alpha}{({{N{({{\mathbb{e}}^{- {\alpha\tau}} - {\mathbb{e}}^{- {\alpha\upsilon}}})}} - {{({1 - {\mathbb{e}}^{- {\alpha\delta}}})}{\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}{\mathbb{e}}^{{- \alpha}\quad t_{n,k}}}}}})}}}}} = {A^{K + m_{A} - 1}\alpha^{m_{\alpha} - 1}e^{- {A{({r_{A} + {\frac{1}{\alpha}{({{N{({{\mathbb{e}}^{- {\alpha\tau}} - {\mathbb{e}}^{- {\alpha\upsilon}}})}} - {{({1 - {\mathbb{e}}^{- {\alpha\delta}}})}{\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}{\mathbb{e}}^{{- \alpha}\quad t_{n,k}}}}}})}}})}}}e^{- {\alpha{({r_{\alpha} + {\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}t_{n,k}}}})}}}}}} & {{Equation}\quad 11}\end{matrix}$and hence $\begin{matrix}{\begin{matrix}{{P\left( a \middle| D \right)} = {\int{{P\left( {A,\left. a \middle| D \right.} \right)}{dA}}}} \\{\propto {\int{A^{K + m_{A} - 1}a^{m_{a} - 1}e^{- {A{({r_{A} + {\frac{1}{a}{({{N{({e^{- {a\tau}} - e^{- {av}}})}} - {{({1 - e^{- {a\delta}}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{- {at}_{n,k}}}}}})}}})}}}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}{dA}}}} \\{= {a^{m_{a} - 1}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}{\int{A^{K + m_{A} - 1}e^{- {A{({r_{A} + {\frac{1}{a}{({{N{({e^{- {a\tau}} - e^{- {av}}})}} - {{({1 - e^{- {a\delta}}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{- {at}_{n,k}}}}}})}}})}}}{dA}}}}} \\{= \frac{a^{m_{a} - 1}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}{\Gamma\left( {K + m_{A}} \right)}}{\left( {r_{A} + {\frac{1}{a}\left( {{N\left( {e^{- {a\tau}} - e^{- {av}}} \right)} - {\left( {1 - e^{- {a\delta}}} \right){\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{- {at}_{n,k}}}}}} \right)}} \right)^{K + m_{A}}}}\end{matrix}{while}} & {{Equation}\quad 12} \\{{{P\left( {\left. A \middle| a \right.,D} \right)} = {\frac{q^{K + m_{A}}}{\Gamma\left( {K + m_{A}} \right)}A^{K + m_{A} - 1}e^{- {Aq}}}}{where}{q = {r_{A} + {\frac{1}{a}\left( {{N\left( {e^{- {a\tau}} - e^{- {av}}} \right)} - {\left( {1 - e^{- {a\delta}}} \right){\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{- {at}_{n,k}}}}}} \right)}}}} & {{Equation}\quad 13}\end{matrix}$so that again the conditional posterior distribution for A given a is aGamma distribution.

With a slight modification algorithm A1 can still be used:

In step 510: 2 P(A,α|D) will instead of equation 8 be given by equation11 taking into account the various deadtimes:

-   520: Calculate    ${P\left( {A,\left. a \middle| D \right.} \right)} \propto {A^{K + m_{A} - 1}a^{m_{a} - 1}e^{- {A{({r_{A} + {\frac{1}{a}{({{N{({e^{- {a\tau}} - e^{- {av}}})}} - {{({1 - e^{- {a\delta}}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{- {at}_{n,k}}}}}})}}})}}}e^{- {a{({r_{a} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad t_{n,k}}}})}}}}$    Algorithm A2, Multiple Fluorophores, Unknown Fluorophore Lifetimes,    with Deadtimes.

Now we repeat the analysis, introducing the remaining complication thatthere are several fluorophores, with parameters given by A_(i), α_(i)for i=1, . . . ,I. The deadtimes are independent of which fluorophoresthe photon came from.

We then have${{P\left( {\left. D_{n} \middle| A \right.,a} \right)} = {\left( {\prod\limits_{k = 1}^{K_{n}}\quad\left( {\prod\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}} \right)} \right){P\left( {\left. {{no}\quad{other}\quad{photons}} \middle| A \right.,a} \right)}}},$where A denotes the vector of the A_(i)s, and α denotes the vector ofα_(i)s. As before, we have${P\left( {\left. {{no}\quad{other}\quad{photons}} \middle| A \right.,a} \right)} = {e^{- {({{\int_{r}^{t_{n,1}}{{\lambda{(t)}}\quad{\mathbb{d}t}}} + {\sum\limits_{k = 1}^{K_{n} - 1}\quad{\int_{t_{n,k}}^{t_{n,{k + 1}}}{{\lambda{(t)}}\quad{\mathbb{d}t}}}} + {\int_{t_{n,K_{n}}}^{v}{{\lambda{(t)}}\quad{\mathbb{d}t}}}})}} \approx e^{- {({{\int_{r}^{v}{{\lambda{(t)}}\quad{\mathbb{d}t}}} - {\sum\limits_{k = 1}^{K_{n}}\quad{\int_{t_{n,k}}^{t_{n,k} + \delta}{{\lambda{(t)}}\quad{\mathbb{d}t}}}}})}}}$where the last line follows if we make the assumption that none of thedeadtimes overlap with the period following υ. Thus $\begin{matrix}{{{P\left( {\left. {{no}\quad{other}\quad{photons}} \middle| A \right.,a} \right)} \approx e^{- {({{\int_{r}^{v}{{\lambda{(t)}}\quad{\mathbb{d}t}}} - {\sum\limits_{k = 1}^{K_{n}}\quad{\int_{t_{n,k}}^{t_{n,k} + \delta}{{\lambda{(t)}}\quad{\mathbb{d}t}}}}})}}} = e^{- {({{\int_{r}^{v}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t}\quad{\mathbb{d}t}}}} - {\sum\limits_{k = 1}^{K_{n}}\quad{\int_{t_{n,k}}^{t_{n,k} + \delta}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t}\quad{\mathbb{d}t}}}}}})}}} \\{= {e^{\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({{\lbrack e^{{- a_{i}}t}\rbrack}_{r}^{v} - {\sum\limits_{k = 1}^{K_{n}}\quad{\lbrack e^{{- a_{i}}T}\rbrack}_{r_{k,n}}^{v_{n,k} + \delta}}})}}} = e^{- {\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})} - {\sum\limits_{k = 1}^{K_{n}}\quad{({e^{{- a_{i}}t_{n,k}} - e^{- {a_{i}{({t_{n,k} + \delta})}}}})}}})}}}}}} \\{= {e^{- {\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})} - {\sum\limits_{k = 1}^{K_{n}}\quad{(e^{{- a_{i}}{t_{n,k}{({1 - e^{{- a_{i}}\delta}})}}})}}})}}}} = e^{- {\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}})}}}}}}\end{matrix}$so combining all the data from N excitations together (independent fromexcitation to excitation), we have $\begin{matrix}{{P\left( {\left. D \middle| A \right.,a} \right)} = {\prod\limits_{n = 1}^{N}\quad\left( {\left( {\prod\limits_{k = 1}^{K_{n}}\quad\left( {\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}} \right)} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}})}}}}} \right)}} \\{= {\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}})}}}}}}} \\{= {\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{{{- N}{\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}}}} + {\sum\limits_{i = 1}^{I}\quad{({\frac{A_{i}}{a_{i}}{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}})}}}}}\end{matrix}.$

Now, applying Bayes' theorem to determine the posterior, we find that$\begin{matrix}\begin{matrix}{{P\left( {A,\left. a \middle| D \right.} \right)} = \frac{{P\left( {A,a} \right)}{P\left( {\left. D \middle| A \right.,a} \right)}}{\int{{P\left( {A,a} \right)}{P\left( {\left. D \middle| A \right.,a} \right)}{dAda}}}} \\{\propto {\prod\limits_{i = 1}^{I}\quad{\frac{r_{A_{i}}^{m_{A_{i}}}A_{i}^{m_{A_{i}} - 1}e^{{- r_{A_{i}}}A_{i}}}{\Gamma\left( m_{A_{i}} \right)}{\prod\limits_{i = 1}^{I}\quad\frac{r_{a_{i}}^{m_{a_{i}}}a_{i}^{m_{a_{i}} - 1}e^{{- r_{a_{i}}}a_{i}}}{\Gamma\left( m_{a_{i}} \right)}}}}} \\{\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{{{- N}{\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}}}} + {\sum\limits_{i = 1}^{I}\quad{({\frac{A_{i}}{a_{i}}{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}})}}}} \\{\propto {\left( {\prod\limits_{i = 1}^{I}\quad{A_{i}^{m_{A_{i}} - 1}a_{i}^{m_{a_{i}} - 1}}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)}} \\{e^{{- {\sum\limits_{i = 1}^{I}\quad{r_{A_{i}}A_{i}}}} - {\sum\limits_{i = 1}^{I}\quad{r_{a_{i}}a_{i}}} - {N{\sum\limits_{i = 1}^{I}\quad{\frac{A_{i}}{a_{i}}{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}}}} + {\sum\limits_{i = 1}^{I}\quad{({\frac{A_{i}}{a_{i}}{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}})}}}} \\{= {\left( {\prod\limits_{i = 1}^{I}\quad{A_{i}^{m_{A_{i}} - 1}a_{i}^{m_{a_{i}} - 1}}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {r_{a_{i}}a_{i}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}}}\end{matrix} & {{Equation}\quad 14}\end{matrix}$

We now use algorithm A2, which though it appears similar to A1 differsin step 520 where P(A,α|D) will be given, instead of by equation 8, byequation 14 taking into account the various deadtimes and thepossibility of a plurality of unknown fluorophore lifetimes and/orintensities.

-   510: 2 Calculate, for each point on the multidimensional grid,    ${P\left( {A,\left. a \middle| D \right.} \right)} = {\left( {\prod\limits_{i = 1}^{I}\quad{A_{i}^{m_{A_{i}} - 1}a_{i}^{m_{a_{i}} - 1}}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {r_{a_{i}}a_{i}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}}$

As can be understood both from the appearance of equation 14 and, moreimportantly, from the fact that the grid is now multidimensional, thisalternative of algorithm A1 will often be computational heavy andtime-consuming.

Algorithm A3, Multiple Fluorophores with Known Lifetimes

A third embodiment of the inventive method of the present inventionrefers to measurement cases in which the lifetimes of the fluorophoresin the sample are known. This can be the case for example when complexorganic molecules are “tagged” with known fluorophores. The wantedmeasure may then be for example the quantities of the fluorophores,reflected in the initial intensities A_(i), or the ratios between thequantities of the different fluorophores, reflected in the ratios of theinitial intensities A_(i)

Equation 14 above is simplified in the case of known lifetimes:${P\left( A \middle| D \right)} \propto {\left( {\prod\limits_{i = 1}^{I}\quad A_{i}^{m_{A_{i}} - 1}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}}$but this is not sufficient to get a simple fast implementation, becauseof the number of terms in the second factor. By, in addition, inferringwhich fluorophore each photon came from it will be shown that thecomputation can be simplified, and hence the time needed for theBayesian inference analysis is significantly shorter. Surprisingly, wegain here by introducing a large number of extra variables—in fact weintroduce I new variables per photon received. For each photon φ_(n,k)received we introduce the I variables s_(n,k,i), which is defined to be1 if φ_(n,k) is emitted by fluorophore i and is 0 otherwise. The svariable tells which fluorophore a photon originates from. Let S denotethe Cartesian product of all these new variables.

We can then recomplicate the above equation as follows: $\begin{matrix}{{P\left( A \middle| D \right)} = {{\int{{P\left( {A,\left. S \middle| D \right.} \right)}{dS}}} = \frac{\int{{P(A)}{P\left( D \middle| A \right)}{P\left( {\left. S \middle| D \right.,A} \right)}{dS}}}{P(D)}}} \\{\propto {\int{\left( {\prod\limits_{i = 1}^{I}\quad A_{i}^{m_{A_{i}} -}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}{P\left( {\left. S \middle| D \right.,A} \right)}{dS}}}}\end{matrix}$ ${Now},\begin{matrix}{{P\left( {\left. S \middle| D \right.,A} \right)} = {{\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{P\left( {\left. s_{n,k} \middle| D \right.,A} \right)}}} = {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{P\left( {\left. s_{n,k} \middle| {{photon}\quad{referred}\quad{to}\quad{by}\quad s_{n,k}\quad{is}\quad{at}\quad t_{n,k}} \right.,A} \right)}}}}} \\{= {{\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad\frac{P\left( {{there}\quad{is}\quad a\quad{photon}\quad{from}\quad{fluor}\quad i\quad{at}\quad t_{n,k}} \middle| A \right)}{P\left( {{there}\quad{is}\quad a\quad{photon}\quad{from}\quad{fluor}\quad i\quad{at}\quad t_{n,k}} \middle| A \right)}}} = {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad\frac{\sum\limits_{i = 1}^{I}\quad{S_{n,k,i}A_{i}e^{{- a_{i}}t_{n,k}}}}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}}}}\end{matrix}$ therefore $P\begin{matrix}{\left( A \middle| D \right) \propto {\int{\left( {\prod\limits_{i = 1}^{I}\quad A_{i}^{m_{A_{i}} - 1}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad\left( {\frac{\sum\limits_{i = 1}^{I}\quad{s_{n,k,i}A_{i}e^{{- a_{i}}t_{n,k}}}}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}} \right)}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}{dS}}}} \\{= {\int{\left( {\prod\limits_{i = 1}^{I}\quad A_{i}^{m_{A_{i}} - 1}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{s_{n,k,i}A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}{dS}}}}\end{matrix}$

By simply removing the integration throughout the above argument, we canalso arrive at $\begin{matrix}{{P\left( {A,\left. S \middle| D \right.} \right)} \propto {\left( {\prod\limits_{i = 1}^{I}\quad A_{i}^{m_{A_{i}} - 1}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{s_{n,k,i}A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}}} & {{Equation}\quad 15}\end{matrix}$which is the starting point for algorithm A3.

We now plan to take random samples from the joint distribution of A andS. In order to do so we will use Markov Chain Monte Carlo techniques. Inparticular, given a random sample (A_(j),S_(j)) we will, if j is even,set S_(j+1) equal to S_(j), and draw A_(j+1), from P(A|S_(j+1)), and ifj is odd, set A_(j+1) equal to A_(j) and draw S_(J+1), fromP(S|A_(j+1)).

Drawing S_(j+1) from P(S↑A_(j+1)) is easy, because this conditionaldistribution of S is separable over its K components S_(n,k), and eachof these is a simple discrete distribution: $\begin{matrix}\begin{matrix}{{P\left( {\left. S \middle| A \right.,D} \right)} \propto {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{s_{n,k,i}A_{i}e^{{- a_{i}}t_{n,k}}}}}}} \\{{\therefore{P\left( {{s_{n,k,i} = 1},{{\left( {\forall{i^{\prime} \neq i}} \right)s_{n,k,i^{\prime}}} = \left. 0 \middle| A \right.},D} \right)}} = \frac{A_{i}e^{{- a_{i}}t_{n,k}}}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}\end{matrix} & {{Equation}\quad 16}\end{matrix}$and each s_(n,k) may be handled independently.

Drawing A_(j+1, from P(A|S) _(j+1)) requires sampling from 1 Gammadistributions; how to do this will be exemplified in an appendix. Thedistributions in question are given by: $\begin{matrix}\begin{matrix}{{P\left( {\left. A \middle| S \right.,D} \right)} \propto {\left( {\prod\limits_{i = 1}^{I}\quad A_{i}^{m_{A_{i}} - 1}} \right)\left( {\prod\limits_{i = 1}^{I}\quad{\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}}} \\{{\therefore{P\left( {\left. A_{i} \middle| S \right.,D} \right)}} = {A_{i}^{m_{A_{i}} - 1 + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad s_{n,k,i}}}}e^{- {A_{i}{({r_{A_{i}} + {\frac{1}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}}}\end{matrix} & {{Equation}\quad 17}\end{matrix}$

Having got our sequence of random samples (A₁,S₁(A₂,S₂₁ . . . we convertit into a sequence of samples of A by discarding the values of S.

With the theoretical background above, and with reference to theflowchart of FIG. 3, an algorithm suitable for software implementationwill be described. The algorithm will be denoted A3.

Algorithm A3:

-   600: In a first initialising step 600, which is included in step 300    of the flowchart of FIG. 3, the user enters functions, values or    distributions describing a and A.-   600: 1 a) If based on information from previous experiments:-   600: 1.1 i) and if posterior of previous experiments are of Gamma    form (eq. 3):    -   copy m and r from the posterior of the previous experiment-   600: 1.2 ii) if other distribution form, e.g. f(A_(i)):    -   plot f(A_(i)) and adjust m_(A) _(i) ,r_(A) _(i) until Gamma        functions mimic f(A_(i)) or    -   consider information loss        $I_{m,r} = {\int{{f(A)}\log\frac{f(A)}{\Gamma_{m,r}(A)}{dA}}}$        and adjust m and r    -   for minimum I_(m,r) using some minimisation routine.-   600: 2 b) If information is vague knowledge of fluorophore; for each    fluorophore pick m_(A) _(i) ,r_(A) _(i) , so that distribution fits    “mental image”:-   600: 2.1 i) for each A_(i) do-   600: 2.1.1 a) pick m_(A) _(i) >0 (that the larger the value the    narrower the prior)-   600: 2.1.2 b) pick mean μ of prior (e.g. 10⁸ meaning 1 photon per 10    ns)-   600: 2.1.3 c) set r_(A) _(i) =m_(A) _(i) /μ-   600: 2.1.4 d) plot prior on screen (linear or log scale)-   600: 2.1.5 e) Does the plot fit the mental image desired ? Yes—new    A_(i)(goto i), No—repeat (a) to (e) as necessary-   605: In a measuring step 605, which corresponds to the steps    305-320, datasets of arrival times are collected D_(n).-   610: The Bayesian inference analysis, step 610 (corresponding    to 325) comprises the substeps (610: 1-4) of:-   610: 1 Initialize:-   610: 1.1 set number of samples required-   610: 1.2 select a starting value for each 4, options:    -   i) select random sample from prior on A_(i)    -   ii) select mean from prior on A_(i)    -   iii) select mode from prior on A_(i)    -   iv) user input data-   610: 2 Update S(s_(n,k) representing the source of each photon)-   610: 2.1 For each n(1,2, . . . ,K_(n))-   610: 2.1.1 For each k(1,2, . . . , K_(n))-   610: 2.1.1.1 Update s_(n,k):    -   select a random sample from the discrete distribution (eq. 17):        ${P\left( {{s_{n,k,i} = 1},{{\left( {\forall{i^{\prime} \neq i}} \right)s_{n,k,i^{\prime}}} = \left. 0 \middle| A \right.},D} \right)} = \frac{A_{i}e^{{- a_{i}}t_{n,k}}}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}$-   610: 2.1.2 Next k-   610: 2.2 Next n-   610: 3 Update A-   610: 3.1 For each i(1,2, . . . , I)-   610: 3.1.1 Update A_(i):-   610: 3.1.1.1 Calculate    $m = {m_{A_{i}} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad s_{n,k,i}}}}$-   610: 3.1.1.2 Calculate    $r_{i} = {r_{A_{i}} + {\frac{1}{a_{i}}\left( {{N\left( {e^{{- a_{i}}\tau} - e^{{- a_{i}}v}} \right)} - {\left( {1 - e^{{- a_{i}}\delta}} \right){\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}} \right)}}$-   610: 3.1.1.3 Set new value of A_(i) to be a random sample from the    Gamma distribution with parameters m, r_(i) (see Appendix)-   610: 3.1 Next i-   610: 4 Add value of A to list of samples of A-   610: 5 Enough samples? Yes—goto 620, No—repeat 610: 2-5-   620: In an outputting step 620 (corresponding to 330) the resulting    list of samples is presented to the user. Optionally histograms    (corresponding to the marginal distributions of previous algorithms)    can be presented. Other measures can be calculated from the    histograms for example the marginal means with the 5^(th) and    95^(th) centiles.

Alternatively it is equally feasible to start with the source variablesS. The algorithm would be slightly modified, the important modificationsare:

-   610: 1 Initialize:    -   select some value at random for each s′_(n,k) defined by        $s_{n,k}^{\prime} \equiv {\sum\limits_{i = 1}^{I}\quad{is}_{n,k,i}}$        from the set {1,2, . . . ,I}        and-   610: 2 Update A-   610: 3 Update S

Note that the procedure, in step 610: 3, to draw a random sample fromthe Gamma distribution with parameters m,r_(i) can be done in variousways. A suitable method will be described in the Appendix.

FIG. 5 illustrates three prior distributions entered by a user insubstep 610:1. In FIG. 6 is a typical output from step 620 shown. FIG. 6shows histograms for three different fluorophores. The x-marked linesindicate true values.

As apparent to the skilled in the art the above algorithm can beconverted to software code in many different ways and care should betaken to achieve computationally effective and efficient code, e.g. inthe loops calculate only the parts that change in the loop. Such, andother simplifications, should be obvious to the skilled programmer andnot considered to be a part of the invention.

Algorithms A4, A5, and A6: Multiple Fluorophores with Known Lifetimesand Binning of the Photon Times

In a further embodiment of the inventive method of the present inventionthe previous algorithm (A3) is modified to also account for the factthat the measurement system has a finite resolution in time. Theresolution of the instrument is typically given by the sampling rate ofthe photon timing unit 150, which could be of the order of 2G samples/s,giving a resolution of 0.5 ns. The finite resolution has the effect thatit is not possible to separate in time photons arriving less than thetime resolution apart; they are “binned” together. Instead of a set ofarrival times, the data really consist of a fixed set of bin positionsand for each excitation the number of photons arriving in each bin. Thisbinning in time of the arrival times of the photons can beadvantageously utilized in the Bayesian inference analysis; indeedalgorithm A4 may in some ways be more useful than A3 even when the dataare not binned, as it is easy to bin the data before using an inferencealgorithm.

Considering the procedure of algorithm A3, and in particular the stepwhen we resample A_(j+1) from P(A|S_(j+1)), we note that the onlydynamic information we need to do this, other than the photon arrivaltimes, is the set of quantities${\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}s_{n,k,i}}},$as we can see from the resampling equation (which defines the parametersof a Gamma distribution): $\begin{matrix}{{P\left( {\left. A_{i} \middle| S \right.,D} \right)} = {A_{i}^{m_{A_{i}} - 1 + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad s_{n,k,i}}}}{e^{- {A_{i}{({r_{A_{i}} + {\frac{1}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}.}}} & {{Equation}\quad 18}\end{matrix}$

However, this set of quantities simply tells us how many photons arebelieved, by the current posterior sample, to come from eachfluorophore.

Similarly, when resampling S_(j+1) from P(S|A_(j+1)), the only dynamicinformation we need is the value of A_(j+1), as can be seen from theresampling equation (which defines a discrete distribution):$\begin{matrix}{{P\left( {{s_{n,k,i} = 1},{{\left( {\forall{i^{\prime} \neq i}} \right)s_{n,k,j^{\prime}}} = \left. 0 \middle| A \right.},D} \right)} = {\frac{A_{i}e^{{- a_{i}}t_{n,k}}}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}.}} & {{Equation}\quad 19}\end{matrix}$

Thus it doesn't matter whether we remember or forget which photon iscurrently believed to come from which fluorophore; simply rememberinghow many there are in each time bin, and how many (total) are allocatedto each fluorophore, is sufficient.

Instead of keeping track of, and resampling, the S variables then, weneed only remember, and resample, the number of photons in each time binbelieved to come from each fluorophore.

Formally then, we define b_(k,i) to be the number of photons in time bink that came from fluorophore i, taking all the excitations together. Letc_(n,k) be the total number of photons in time bin k of excitation n;then C constitutes the data. Let d_(k) be the midtime of bin k and w_(k)be its width; we assume that for all other t in that same bine^(−max_(i)(a_(i))_(t))differs negligibly from $e^{- {\max\limits_{i}{{(a_{i})}d_{k}}}}$(and if this is not so we will use algorithm A7 below instead).

As can be seen from the resampling equation for S_(j+1) therefore, theconditional distribution of b_(k) given A_(j+1) will be multinomial,with multi-way toss probabilities given by${{{EMBED}\frac{w_{k}A_{i}e^{{- a_{i}}d_{k}}}{\sum\limits_{i = 1}^{I}\quad{w_{k}A_{i}e^{{- a_{i}}d_{k}}}}} = \frac{A_{i}e^{{- a_{i}}d_{k}}}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}d_{k}}}}},$and total number of tosses $\begin{matrix}{\sum\limits_{n = 1}^{N}\quad{c_{n,k}.}} & {{Equation}\quad 20}\end{matrix}$

For odd j we now draw a sample from a multinomial distribution with themulti-way toss probabilities defined above, while for even j we resampleA_(j+1) according to $\begin{matrix}{{P\left( {\left. A_{i} \middle| B \right.,D} \right)} = {A_{i}^{m_{A_{i}} - 1 + {\sum\limits_{k = 1}^{K}\quad b_{k,i}}}e^{- {A_{i}{({r_{A_{i}} + {\frac{1}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K}\quad{c_{n,k}e^{{- a_{i}}d_{k}}}}}}})}}})}}}}} & {{Equation}\quad 21}\end{matrix}$

How to draw samples from a multinomial distribution efficiently will beshown in appendix. Asymptotically it takes only as long as the squareroot of the total number of photons in the bin multiplied by the largesttoss probability. The run time is no longer, as in algorithm A3,proportional to the number of photons received, but to the square rootof that number.

With the above reasoning and derivation algorithm A3 can be modifiedaccordingly and is described below and with references to the flowchartsof FIG. 3.

Algorithm A4:

-   700: In a first initialising step 700, which is included in step 300    of the flowchart of FIG. 3, the user enters functions, values or    distributions describing α and A.-   700: 1 a) If information from previous experiments:-   700: 1.1 i) and if posterior of previous experiments are of Gamma    form (eq. 3):    -   copy m and r from the posterior of the previous experiment-   700: 1.2 ii) if other distribution form, e.g. f(A_(i)):    -   plot f(A_(i)) and adjust m_(A) _(i) ,r_(A) _(i) until Gamma        functions mimics f(A_(i)) or    -   consider information loss        $I_{m,r} = {\int{{f(A)}\log\frac{f(A)}{\Gamma_{m,r}(A)}{dA}}}$        and adjust m and r for minimum I_(m,r) using some minimisation        routine.-   700: 2 b) If information is vague knowledge of fluorophore; for each    fluorophore pick m_(A) _(i) ,r_(A) _(i) r so that distribution fits    “mental image”:-   700: 2.1 i) for each A_(i) do-   700: 2.1.1 a) pick m_(A) _(i) >0 (that the larger the value the    narrower the prior)-   700: 2.1.2 b) pick mean μ of prior (e.g. 1080-1 photon per 10 ns)-   700: 2.1.3 c) set r_(A) _(i) =m_(A) _(i) /μ-   700: 2.1.4 d) plot prior on screen (linear or log scale)-   700: 2.1.5 e) does the plot match the desired mental image ? Yes—new    A_(i) (goto i), No—repeat (a) to (e) as necessary-   705: In a measuring step 705, which corresponds to the steps    305-320, datasets of arrival times are collected D_(n).-   710: The Bayesian inference analysis, step 710 (corresponding    to 325) comprises the substeps (710: 1-4) of:-   710: 1 Initialize:-   710: 1.1 set bin size and position (from hardware design)-   710: 1.2 count photons in each time bin:    $c_{k} = {\sum\limits_{n = 1}^{N}\quad c_{n,k}}$-   710: 1.3 set number of samples required-   710: 1.4 select starting value for each A_(i), options:    -   i) select random sample from prior on A_(i)    -   ii) select mean from prior on A_(i)    -   iii) select mode from prior on A_(i)    -   iv) user input data-   710: 2 Update B-   710: 2.1 For each k(1,2, . . . ,K)-   710: 2.1.1 Update b_(k)-   710: 2.1.1.1 set b_(k) to be a random sample from the multinomial    distribution with $c_{k} = {\sum\limits_{n = 1}^{N}\quad c_{n,k}}$    total tosses and toss probabilities    ${Pi} = {\frac{A_{i}e^{{- a_{i}}d_{k}}}{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}d_{k}}}}\left( {{see}\quad{Appendix}} \right)}$-   710: 2.2 Next k-   710: 3 Update A-   710: 3.1 For each i(1,2, . . . ,I)-   710: 3.1.1 Update A_(i)-   710: 3.1.1.1 Calculate    $m = {m_{A_{i}} + {\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K}\quad b_{k,i}}}}$-   710: 3.1.1.2 Calculate    $r_{i} = {r_{A_{i}} + {\frac{1}{a_{i}}\left( {{N\left( {e^{{- a_{i}}\tau} - e^{{- a_{i}}v}} \right)} - {\left( {1 - e^{{- a_{i}}\delta}} \right){\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K}\quad{c_{k}e^{{- a_{i}}d_{k}}}}}}} \right)}}$-   710: 3.1.1.3 Set new value of A_(i) to be a random sample from the    Gamma distribution with parameters m, r_(i) (see Appendix)-   710: 3.2 Next i-   710: 4 Add value of A to list of samples of A-   710: 5 Enough samples? Yes—goto 720, No—repeat 710: 1-5-   720: In an outputting step 720 (corresponding to 330) the resulting    list of samples is presented to the user. Optionally histograms    (corresponding to the marginal distributions of previous algorithms)    can be presented. Other measures can be calculated from the    histograms for example the marginal means with the 5^(th) and    95^(th) centiles and/or fractions of total intensity.

Similarly to previous algorithm it is possible to reverse the order ofupdating A and B with slight modifications of the algorithm.

The majority of the running time of algorithm A4 is occupied in drawingrandom samples from multinomial distributions (substep 710: 2).Replacing those samples with their mean will result in a fasteralgorithm. A noticeable deterioration in accuracy of the result can beexpected, although under some circumstances it may still be consideredusable. The algorithm modified in this way will be referred to as A5, Afurther simplification is to replace the Gamma samples (substep 710: 3)with their means, which will result in an even faster but veryinaccurate algorithm, which will be referred to as A6.

Algorithm A7. Multiple Fluorophore with Known Lifetimes and Binning ofthe Photon Times in Bins and Subbins

A further embodiment of the inventive method of the present inventionillustrates the power and usefulness of the method. Suppose the binsintroduced in the previous algorithm are too large, i.e. the resolutionof the instrument is lower than apparently needed for the specificmeasurement. The algorithm described below, will by introducing subbins,dividing each bin into a number of subbins, make it possible to extractinformation that would otherwise been obscured by the resolution of theinstrument. The subbins do not necessarily have to be of equal size, forexample they could be spaced logarithmically; nor is it necessary thateach bin contains the same number of subbins as each other bin.

We now introduce a number of extra variables similar to the B variables,but indicating not only which fluorophore each photon comes from, butalso which subbin it falls in within the bin in which it is known tofall.

Resanpling of A_(j+1) from P(A|B_(j,+1)) occurs in exactly the same wayas in A4. Resampling of B_(j+1) from P(B|A_(j+1)), however, now occursover, for each bin, a multinomial distribution with the number ofoptions being the product of the number of fluorophores and the numberof subbins, and the total number of tosses being the total number ofphotons in that bin, and the multi-way toss probabilities beingdetermined by the subbin centre time and fluorophore decay time that arerelevant. Finally, of course, we discard the subbin information just aswe discard the B variable information.

In detail then, suppose the original time bins are centered on timesd_(k) as before, but that bin k is now divided into L_(k) subbins withtimes centered on d_(k,l) and of width w_(k,j). Similarly, let b_(k,i,l)denote the number of photons currently supposed to have come fromfluorophore i and to be in subbin l of bin k. Then the multinomialdistribution over which we have to resample b_(k) has now IL_(k)possible toss outcomes, still with $\begin{matrix}{{EMBED},\frac{w_{k,l}A_{i}e^{{- a_{i}}d_{k,j}}e^{\sum\limits_{i^{\prime} = 1}^{I}\quad{({\frac{A_{i^{\prime}}}{a_{i^{\prime}}}{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K}\quad e^{{- a_{i^{\prime}}}d_{k,l}}}})}}}{\sum\limits_{i = 1}^{i}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{w_{k,l}A_{i}e^{{- a_{i}}d_{k,l}}e^{\sum\limits_{i^{\prime} = 1}^{I}\quad{({\frac{A_{i^{\prime}}}{a_{i^{\prime}}}{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K}\quad e^{{- a_{i^{\prime}}}d_{k,l}}}})}}}}}} & {{Equation}\quad 22}\end{matrix}$total tosses. However, calculating the toss probabilities is slightlymore difficult than it was before. Working from equation 20 forP(A,S|D), we find that the probability for toss outcome (i,l) in bin kis $\sum\limits_{n = 1}^{N}\quad c_{n,k}$while when resampling A_(j+1) we must work from $\begin{matrix}{{{P\left( {\left. A_{i} \middle| B \right.,D} \right)} = {A_{i}^{m_{A_{i}} - 1 + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,j}}}}e^{- {A_{i}{({r_{A_{i}} + {\frac{1}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{({e^{{- a_{i}}d_{k,l}}{\sum\limits_{i^{\prime} = 1}^{I}\quad b_{k,l^{\prime},l}}})}}}}})}}})}}}}},} & {{Equation}\quad 23}\end{matrix}$which as usual is a Gamma distribution.

With the above reasoning and derivation, algorithm A4 can be modifiedaccordingly and is described below and with references to the flowchartsof FIG. 3. The modifications occur in the substeps of the Bayesianinference analysis, substeps 810; 1-5 corresponding to step 710: 1-5 ofalgorithm A4, only. Therefore only these substeps are described below.

Algorithm A7:

-   -   810: 1 The Bayesian inference analysis, step 810 (corresponding        to 710) comprises the substeps of:    -   810: 1 Initialize:    -   810: 1.1 set bin size and position: d_(k)(k=1,2, . . . ,K) (from        hardware design)    -   810: 1.2 set subbin size and position for each k: d_(k,j)(l=1,2,        . . . ,L_(k)) (inside d_(k)'s bin)    -   810: 1.3 count photons in each time bin:        $c_{k} = {\sum\limits_{n = 1}^{N}\quad c_{n,k}}$    -   810: 1.4 set number of samples required    -   810: 1.5 select starting value for each A_(i), options:        -   i) select random sample from prior on A_(i)        -   ii) select mean from prior on A_(i)        -   iii) select mode from prior on A_(i)        -   iv) user input data    -   810: 2 Update B    -   810: 2.1 For each k(1,2, . . . , K)    -   810: 2.1.1 Update b_(k):    -   set b_(k) to be a random sample from the multinomial        distribution with        $c_{k} = {\sum\limits_{n = 1}^{N}\quad c_{n,k}}$        total tosses, IL_(k) outcomes and toss probability        ${EMBED}\frac{w_{k,l}A_{i}e^{{- a_{i}}d_{k,j}}e^{\sum\limits_{i^{\prime} = 1}^{I}\quad{({\frac{A_{i^{\prime}}}{a_{i^{\prime}}}{({1 - e^{{- a_{i^{\prime}}}\delta}})}{\sum\limits_{k = 1}^{K}\quad e^{- a_{i^{\prime}d_{k,l}}}}})}}}{{\sum\limits_{i = 1}^{i}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{w_{k,l}A_{i}e^{{- a_{i}}d_{k,l}}}}}e^{\sum\limits_{i^{\prime} = 1}^{I}\quad{({\frac{A_{i^{\prime}}}{a_{i^{\prime}}}{({1 - e^{{- a_{i^{\prime}}}\delta}})}{\sum\limits_{k = 1}^{K}\quad e^{{- a_{i^{\prime}}}d_{k,l}}}})}}}\left( {{see}\quad{Appendix}} \right)$    -   810: 2.2 Next k    -   810: 3 Update A    -   810: 3.1 For each i(1,2, . . . ,I)    -   810: 3.1.1 Update A_(i):

-   810: 3.1.1.1 Calculate    $m = {m_{A_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,l}}}}$    -   810: 3.1.1.2 Calculate        $r_{i} = {r_{A_{i}} + {\frac{1}{a_{i}}\left( {{N\left( {e^{{- a_{i}}\tau} - e^{{- a_{i}}v}} \right)} - {\left( {1 - e^{{- a_{i}}\delta}} \right){\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{\left( e^{{- a_{i}}d_{k}} \right){\sum\limits_{i^{\prime} = 1}^{I}\quad b_{k,i^{\prime},l}}}}}}} \right)}}$    -   810: 3.1.1.3 Set new value of A_(i) to be a random sample from        the Gamma distribution with parameters m,r_(i)(see Appendix)    -   810: 3.2 Next i    -   810: 4 Add value of A to list of samples of A    -   810: 5 Enough samples? Yes—goto 820, No—repeat 810: 2-5

Just as in the previous algorithm it is possible to reverse the order ofupdating A and B with slight modifications of the algorithm.

Algorithm A8, Multiple Fluorophore with Unknown Lifetimes and Binning ofthe Photon Times in Bins and Subbins

The principle of binning of the photon arrival times in bins and subbinscan be utilized also in the case of multiple fluorophores with unknownlifetimes, which will be illustrated below in an embodiment of thepresent invention. It can be seen as a generalization of algorithm A7.As before when the lifetimes are unknown, we presume Gamma priors onthem. The simplifications that result in the simpler case where(∀k)L_(k)=1 should be obvious.

As enunciated in the description of algorithm A2 above, the posteriorhere before binning is given by $\begin{matrix}{{P\left( {A,\left. a \middle| D \right.} \right)} \propto {\left( {\prod\limits_{i = 1}^{I}\quad{A_{i}^{m_{A_{i}} - 1}a_{i}^{m_{a_{i}} - 1}}} \right)\left( {\prod\limits_{n = 1}^{N}\quad{\prod\limits_{k = 1}^{K_{n}}\quad{\sum\limits_{i = 1}^{I}\quad{A_{i}e^{{- a_{i}}t_{n,k}}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {r_{a_{i}}a_{i}} + {\frac{A_{i}}{a_{i}}{({{N{(e^{{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{n = 1}^{N}\quad{\sum\limits_{k = 1}^{K_{n}}\quad e^{{- a_{i}}t_{n,k}}}}}})}}})}}}}} & {{Equation}\quad 24}\end{matrix}$while with binned data and the additional variables, using the notationof the last algorithms, we have $\begin{matrix}{{{EMBED}\left( {A,a,\left. B \middle| D \right.} \right)} \propto {\left( {\prod\limits_{k = 1}^{K}\quad{\prod\limits_{l = 1}^{L_{k}}\quad w_{k,l}^{\sum\limits_{i = 1}^{I}\quad b_{k,i,l}}}} \right){\left( {\prod\limits_{i = 1}^{I}\quad{A_{i}^{m_{A_{i}} - 1 + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,l}}}}a_{i}^{m_{a_{i}} - 1}}} \right) \cdot e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {a_{i}{({r_{a_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{b_{k,i,l}d_{k,l}}}}})}} + {\frac{A_{i}}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{({e^{{- a_{i}}d_{k,l}}{\sum\limits_{i^{\prime} = 1}^{I}\quad b_{k,i^{\prime},l}}})}}}}})}}})}}}}}} & {{Equation}\quad 25}\end{matrix}$under the constraint that for all k we have${\sum\limits_{l = 1}^{L_{k}}\quad{\sum\limits_{i = 1}^{I}\quad b_{k,i,l}}} = {\sum\limits_{n = 1}^{N}\quad c_{n,k}}$(and otherwiseP(A,α,B|D)=0).

It is now easy to write down the conditional distributions:$\begin{matrix}{{P\left( {\left. A \middle| a \right.,B,D} \right)} \propto {\left( {\prod\limits_{i = 1}^{I}\quad A_{i}^{m_{A_{i}} - 1 + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,l}}}}} \right)e^{- {\sum\limits_{i = 1}^{I}\quad{({A_{i}{({r_{A_{i}} + {\frac{1}{a_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K}\quad{({\sum\limits_{l = 1}^{L_{k}}\quad{e^{{- a_{i}}d_{k,l}}{\sum\limits_{i^{\prime} = 1}^{I}\quad b_{k,i^{\prime},l}}}})}}}})}}})}})}}}}} & {{Equation}\quad 26}\end{matrix}$(which is a product of Gamma distributions);

P(B|A,α,D) is a product of multinomials, with P(b_(k)|A,α,D) having tossprobabilities $\begin{matrix}{{EMBED}\frac{w_{k,l}A_{i}e^{{- a_{i}}d_{k,l}}e^{\sum\limits_{i^{\prime} = 1}^{I}\quad{({\frac{A_{i^{\prime}}}{a_{i^{\prime}}}{({1 - e^{{- a_{i^{\prime}}}\delta}})}{\sum\limits_{k = 1}^{K}\quad e^{{- a_{i^{\prime}}}d_{k,l}}}})}}}{\sum\limits_{i = 1}^{i}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{w_{k,l}A_{i}e^{{- a_{i}}d_{k,l}}e^{\sum\limits_{i^{\prime} = 1}^{I}\quad{({\frac{A_{i^{\prime}}}{a_{i^{\prime}}}{({1 - e^{{- a_{i^{\prime}}}\delta}})}{\sum\limits_{k = 1}^{K}\quad e^{{- a_{i^{\prime}}}d_{k,l}}}})}}}}}} & {{Equation}\quad 27}\end{matrix}$and total tosses ${\sum\limits_{n = 1}^{N}\quad c_{n,k}};$while finally $\begin{matrix}\begin{matrix}{{P\left( {A,\left. a \middle| B \right.,D} \right)} \propto \left( {\prod\limits_{i = 1}^{I}\quad{A_{i}^{m_{A_{i}} - 1 + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,l}}}}a_{i}^{m_{a_{i}} - 1}}} \right)} \\{e^{- {\sum\limits_{i = 1}^{I}\quad{({{r_{A_{i}}A_{i}} + {a_{i}{({r_{a_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{b_{k,i,l}d_{k,l}}}}})}} + {\frac{A_{i}}{\sigma_{i}}{({{N{({e^{{- a_{i}}\tau} - e^{{- a_{i}}v}})}} - {{({1 - e^{{- a_{i}}\delta}})}{\sum\limits_{k = 1}^{K}\quad{({\sum\limits_{l = 1}^{L_{k}}\quad{e^{{- a_{i}}d_{k,l}}{\sum\limits_{i^{\prime} = 1}^{I}\quad b_{k,i^{\prime},l}}}})}}}})}}})}}}} \\{\therefore{{P\left( {\left. a \middle| B \right.,D} \right)} \propto {\int{{P\left( {A,\left. a \middle| B \right.,D} \right)}{dA}}}}} \\{= {\prod\limits_{i = 1}^{I}\quad\frac{a_{i}^{m_{a_{i}} - 1}e^{- {a_{i}{({r_{a_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{b_{k,i,l}d_{k,l}}}}})}}}{\Gamma\left( {m_{A_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,l}}}} \right)}}{\left( {r_{A_{i}} + {\frac{1}{a_{i}}\left( {{N\left( {e^{{- a_{i}}\tau} - e^{{- a_{i}}v}} \right)} - {\left( {1 - e^{{- a_{i}}\delta}} \right){\sum\limits_{k = 1}^{K}\quad\left( {\sum\limits_{l = 1}^{L_{k}}\quad{e^{{- a_{i}}d_{k,l}}{\sum\limits_{i^{\prime} = 1}^{I}\quad b_{k,i^{\prime},l}}}} \right)}}} \right)}} \right)^{m_{A_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,l}}}}}}}\end{matrix} & {{Equation}\quad 28}\end{matrix}$where of course the Gamma functions do need to be calculated as they arenot constants.

Given these three conditional distributions, we can now in principleresample as follows. For odd j we now draw a sample for B_(j+1) from amultinomial distribution with the multi-way toss probabilities definedabove, while for even j we resample (A_(j+i), α_(j+1)) (α is a vector)by first sampling a from the marginal conditional P(α_(j+1), |B_(j+1),D), and then sampling A_(j+1) from the conditional distributionP(A_(j+1)|α_(j+1),B_(j+1),D).

The only part of the this procedure that has so far not been fullydefined is the step of resampling from P(α_(j+1)|B_(j+1),D). This willbe expanded below.

With the above reasoning and derivation, algorithm A7 can be modifiedaccordingly and is described below and with references to the flowchartsof FIG. 3. The modifications of importance occur in the substeps of theBayesian inference analysis, substeps 910: 1-4 corresponding to substeps810: 1-4 of algorithm A4, only. Therefore only these substeps aredescribed below.

Algorithm A8:

-   -   910: The Bayesian inference analysis, step 910 (corresponding        to 810) comprises the substeps of:    -   910: 1 Initialize:    -   910: 1.1 set bin size and position: d_(k)(k=1,2, . . . ,K) (from        hardware design)    -   910: 1.2 set subbin size and position for each k: d_(k,l)(l=1,2,        . . . ,L_(k)) (inside d_(k)'s bin)    -   910: 1.3 count photons in each time bin:        $c_{k} = {\sum\limits_{n = 1}^{N}\quad c_{n,k}}$    -   910: 1.4 set number of samples required    -   910: 1.5 set both r_(A) _(i) , m_(A) _(i) and r_(α) _(i) ,m_(α)        _(i) in one of the same ways that r_(A) _(i) ,m_(A) _(i) were        set in algorithm A4.    -   910: 2: For each j:    -   910: 2.1 If j is odd: Update B.    -   910: 2.1.1 For each k(1,2, . . . ,K)    -   910: 2.1.1.1 Update b_(k):        -   set b_(k) to be a random sample from the multinomial            distribution with            $c_{k} = {\sum\limits_{n = 1}^{N}\quad c_{n,k}}$            total tosses, IL_(k) outcomes and toss probabilities            ${EMBED}\frac{w_{k,l}A_{i}e^{{- a_{i}}d_{k,l}}e^{\sum\limits_{i^{\prime} = 1}^{I}\quad{({\frac{A_{i^{\prime}}}{a_{i^{\prime}}}{({1 - e^{{- a_{i^{\prime}}}\delta}})}{\sum\limits_{k = 1}^{K}\quad e^{{- a_{i^{\prime}}}d_{k,l}}}})}}}{{\sum\limits_{i = 1}^{I}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{w_{k,l}A_{i}e^{{- a_{i}}d_{k,l}}}}}e^{\sum\limits_{i^{\prime} = 1}^{I}\quad{({\frac{A_{i^{\prime}}}{a_{i^{\prime}}}{({1 - e^{{- a_{i^{\prime}}}\delta}})}{\sum\limits_{k = 1}^{K}\quad e^{- a_{i^{\prime}d_{k,l}}}}})}}}\left( {{see}\quad{Appendix}} \right)$    -   910: 2.1.2 Next k    -   910: 2.2 If j is even:    -   910: 2.2.1 Resample (A_(j+1), α_(j+1)) by    -   910: 2.2.1.1 Sampling a from P(α_(j+1)|B_(j+1),D) (eq. 28) (see        below for description of how this is done)    -   910: 2.2.1.2 Sampling A_(j+1), from P(A_(j+1)|α_(j+1),B_(j+1),D)        (eq.26)    -   910: 2.2.2 Add A to list of samples of A    -   910: 3 next j:    -   910: 4 Enough samples? Yes—goto 920, No—repeat 910: 2-4

Sampling α from P(α_(j+1)|B_(j+1),D) (in step 935) is not obvious. Belowa suitable method, readily converted into program code, is outlined.

Resampling from the conditional for a can be performed as follows: Weintroduce a Markov Chain Monte Carlo (MCMC) technique whose proposaldistribution is not a conditional distribution, does depend on theprevious sample and does not always accept the new sample.

Suppose then that we have values of B_(j+1), and D, and also theprevious value of the ith component of α namely α_(j,i) (we now have twosubscripts in order to indicate both which sample and which component).Note that the conditional distribution of P(α|B,D) is separable over i,so that we can treat each component separately. Let the currentconditional distribution P(α_(i)|B,D) be denoted f(α_(i)), so that$\begin{matrix}{{f\left( a_{i} \right)} \propto \frac{a_{i}^{m_{a_{i}} - 1}e^{- {a_{i}{({r_{a_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad{b_{k,i,l}d_{k,l}}}}})}}}{\Gamma\left( {m_{A_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,l}}}} \right)}}{\left( {r_{A_{i}} + {\frac{1}{a_{i}}\left( {{N\left( {e^{{- a_{i}}\tau} - e^{{- a_{i}}v}} \right)} - {\left( {1 - e^{{- a_{i}}\delta}} \right){\sum\limits_{k = 1}^{K}\quad\left( {\sum\limits_{l = 1}^{L_{k}}\quad{e^{{- a_{i}}d_{k,l}}{\sum\limits_{i^{\prime} = 1}^{I}\quad b_{b,i^{\prime},l}}}} \right)}}} \right)}} \right)^{m_{A_{i}} + {\sum\limits_{k = 1}^{K}\quad{\sum\limits_{l = 1}^{L_{k}}\quad b_{k,i,l}}}}}} & {{Equation}\quad 29}\end{matrix}$

Then the overall move in ai needs to satisfy detailed balance withrespect to this distribution. It is easy to show that if it does, thenthe alternating sequence of moves described before will also holddetailed balance over every subsequence of three moves, and hence thatthe necessary conditions for our MCMC process will be met.

It is also easy to show that if we have two potential types of move thateach holds detailed balance with respect to this distribution, then soalso will a compound move that proceeds by choosing one of these twotypes at random (with any fixed probabilities) and applying it. We nowthen go on to describe two types of move, which we will use with somefixed probabilities (for example 0.5 and 0.5) whenever a move toresample α_(i) is required.

Both types of move involve sampling a proposed value α′ from a proposaldistribution g(α′|α_(j,i)), dependent on the previous value of α_(i),namely α_(j,i), followed by a probabilistic decision either to accept α′(in which case we set a_(j+1,i)=α′) or to reject it (in which case weset α_(j+1,i)=α_(j,i)). In both cases we accept with probability$\begin{matrix}{p_{a} = {\frac{{f\left( a^{\prime} \right)}{g\left( a_{j,i} \middle| a^{\prime} \right)}}{{f\left( a_{j,i} \right)}{g\left( a^{\prime} \middle| a_{j,i} \right)}}.}} & {{Equation}\quad 30}\end{matrix}$

All that remains then is to define the two proposal distributions, onefor each type of move.

Proposal Distribution 1—“The Coarse Grid”:

Let α_(j,i) be given; fix also (once and for all) some parameters γ andH, where γ is a positive real number (20 will do) and H is a positiveinteger (again, 20 will do). Then let X₁ be the set$\left\{ {x_{h} = {a_{j,i}{\gamma^{\frac{h}{H}}:{{h\varepsilon}\left\{ {{- H},{{- H} + 1},\ldots,{H - 1},H} \right\}}}}} \right\},$and let g₁(α′|α_(j,i)) be the unique (discrete) probability distributionon X, that is proportional to f(x).Proposal Distribution 2—“The Fine Fill-In”:

Using the same parameters γ and H, let X₂ be the closed interval$\left\lbrack {{a_{j,i}\gamma^{- \frac{1}{H}}},{a_{j,i}\gamma^{+ \frac{1}{H}}}} \right\rbrack,$and let g₂ (α′|α_(j,i)) be the unique (continuous) probabilitydistribution on X₂ that is proportional to $\frac{1}{x}.$

The measurement system with the capability of detecting and recordingmultiple photons after one excitation in combination with the method ofthe present invention, the essence of which is implementing Bayesianinference in the algorithms executed during each measurement, representsa significant improvement in measurement speed over the prior arttechniques. The inventive method of, in each measurement, inferring thesource (i.e. the originating fluorophore) of each photon, furthershortens the measurement time (algorithm A3). In Table 1 a comparison isgiven of the measurement times of four different dyes (having differentlifetimes) using the TCSPC and the TCMPC techniques. The same solutionswere used with both instruments. These numbers are only indicative, butshow a faster speed for the TCMPC technique. TABLE 1 Measurement timewith different dyes using TCSPC and TCMPC, respectively. Time to run theexperiment(s) Lifetime of the Dye TCSPC TCMPC  4 ns 493 9.94  8 ns  289.55 14 ns 196 18.79 20 ns 640 13.35

In the described method according to the present invention, thefluorescence decay of a fluorophore bound to a biological molecule hasbeen modelled by a mono-exponential decay or in some cases abi-exponential decay. However, there are situations in which one doesnot expect a limited number of discrete decay times but rather adistribution of decay times. In complex heterogeneous biologicalsystems, multi-exponential decays functions can appear to provide abetter fit of the fluorescence decay data than the assumption of amono-exponential decay. But the assumption of multiple discretecomponents is essentially arbitrary and is often erroneous. For example,for a single tryptophan protein, the distribution of decay time mayreflect the distribution of protein conformation, alternatively it couldbe a protein with many tryptophan residues, so it is not practical toconsider the individual decay times. Moreover, interaction between thefluorophore and the environment can result in complex fluorescence decayprofiles that represent a continuous distribution of lifetimes. A methodaccording to the present invention may also account for such cases.

The intensity decays λ(t) may be analysed in terms of a lifetimedistribution ρ(τ_(l), wherein τ_(l)/1α and α is the decay constant. Thetotal decay law becomes: $\begin{matrix}{{\lambda(t)} = {\int_{\tau = 0}^{\infty}{{\rho\left( \tau_{1} \right)}e^{\frac{t}{\tau_{l}}}\quad{\mathbb{d}\tau_{l}}}}} & {{equation}\quad 31}\end{matrix}$

The fluorescence decay of biological systems may in some casesadvantageously be described by a so-called stretched exponential as thestretched exponential reflects a continuous distribution of fluorescencelifetimes. The stretched exponential is introduced by: $\begin{matrix}{{\lambda(t)} = {A_{0}.e^{- {(\frac{t}{\tau_{1}})}^{\frac{1}{h}}}}} & {{equation}\quad 32}\end{matrix}$where h is called the heterogeneity and is related to the distributionof decay times. In the case of a mixture of both substrate and product,each contributing to the fluorescence decay, a sum of two stretchedexponentials may be used: $\begin{matrix}{{\lambda(t)} = {{A_{l}.e^{- {(\frac{t}{\tau_{l_{1}}})}^{\frac{1}{h_{1}}}}} + {A_{2}.e^{- {(\frac{t}{\tau_{l_{2}}})}^{\frac{1}{h_{2}}}}}}} & {{equation}\quad 33}\end{matrix}$

The stretched exponential model may advantageously be combined with theanalysis using Bayesian inference according to the present inventionwhich has been exemplified in the algorithms A1-A8. For examplealgorithm A1, taking into account the various deadtimes and utilizingthe stretched exponential model may be modified as described below.

The Poisson arrival rate of the photons at time t is denoted by λ(t) asdescribed by equation 32. The equations of algorithm A1 (equation 11-13)should be modified accordingly by using the following expression for thelikelihood of observing the data D_(n): $\begin{matrix}{{P\left( {\left. D_{n} \middle| A_{0} \right.,\tau,h} \right)} = {\left( {\prod\limits_{k = 1}^{K_{n}}\quad\left( {\lambda\left( t_{n,k} \right)} \right)} \right)e^{- {({{\int_{\tau}^{v}{{\lambda{(t_{n,k})}}\quad{\mathbb{d}t}}} - {\sum\limits_{k = 1}^{K_{n}}\quad{\int_{t_{n,k}}^{t_{n,k} + \delta}{{\lambda{(t_{n,k})}}\quad{\mathbb{d}t}}}}})}}}} & {{equation}\quad 34}\end{matrix}$where the detector is dead for a duration δ after each photon isreceived, for a duration τ after the excitation, and at all times longerthan υ after excitation. K_(n) is the number of photons receivedfollowing the n^(th) laser pulse and t_(n,k) is the arrival of thek^(th) photon from laser pulse n. Furthermore, the expression for priorprobability distribution should be modified accordingly to included anyprior knowledge on the values of the heterogeneity parameter h. Also theother algorithms, including the algorithms using binning of photonarrival times (A4-A8), may in the same way be modified to includestretched exponentials.

The inventive method, exemplified with algorithms A1-A8, has beendescribed here in combination with the measurement system comprising thephotodetector 160 and the photon timing unit 150 having the ability todetect and output multiple photon arrival times after each excitationlight pulse from the pulsed laser 120. This is considered to be the bestmode of operation. However, as appreciated by those skilled in the art,parts of the method of the present invention can advantageously be usedalso with prior art apparatuses such as the TCSPC and Multiplexed-TCSPC,providing only one (TCSPC) or a limited (Multiplexed-TCSPC) number ofphotons per excitation pulse.

The inventive method and measurement system has been exemplified herewith fluorophores. As appreciated by those skilled in the art the systemand methods are equally applicable in characterising other luminescentsubstances e.g. phosphorescent materials.

Example of luminescent substances that can advantageously be analysedwith the described method and system comprises fluorescent materials(which usually have short lifetimes, typically of the order of <1 ns-100ns) e.g. Fluorescein, CyDyes, and ethidium bromides; and Phosphorescentmaterial (which usually have longer lifetimes (typically a few ms)) e.g.ruthenium chelates, terbium chelates.

From the invention thus described, it will be obvious that the inventionmay be varied in many ways. Such variations are not to be regarded as adeparture from the scope of the invention, and all such modifications aswould be obvious to one skilled in the art are intended for inclusionwithin the scope of the following claims.

APPENDIX

A method of drawing random samples from a Gamma distribution and from amultinomial distribution, respectively, will be given below. Thedescription will be in the form of guidelines which the skilledprogrammer can easily convert into subroutines suitable to call fromalgorithms A3-A8.

Random Samples from a Gamma Distribution

The Gamma distribution is defined by${\Gamma_{m,r}(x)} = \frac{r^{m}x^{m - 1}e^{- {rx}}}{\Gamma(m)}$as described above. Obtaining a random sample from it can be done inmany ways. The best known currently to the inventors proceedsdifferently depending on whether m is less than 1, equal to 1, orgreater than 1. Whichever is being used, we can initially assume thatr=1 and then divide the resulting sample by r afterwards.The Easy Case: m=1

In this case the distribution reduces to an exponential distribution.The cumulative distribution function is thereforeP(x<y)=1−e ^(−y)and we may therefore choose a uniformly distributed random number u fromthe interval [0, 1] and solve the equation u=1−e^(−y) for y which weoutput as our random sample.The Difficult Cases: m<1 and m>1

For both of these cases we use forms of rejection sampling.

Suppose we want to sample from a probability distribution f(x) and arealready able to sample from a probability distribution g(x). Supposemoreover that there exists a constant K such that for all x we have Kf(x)≦g(x). Then it is easy to see that the following algorithm willresult in a sample from f(x) with probability K and with the remainingprobability will result in no sample:

-   -   Draw a random sample y from g(x).    -   With probability Kf(x)/g(x) let y be the resultant sample, and        with the remaining probability let there be no resultant sample.

If, therefore, we keep applying this algorithm until it does yield asample, then we will end up with a sample from f(x) with probability1−(1−K)^(N) after N repeats. Since this probability approaches 1 as Napproaches infinity, we will eventually end up with a sample withprobability 1.

All that remains now is to specify g(x) for each of the relevant cases,state how to draw samples from g(x), and how to calculate Kf (x)/g(x).

m>1

g(x) is given by${g(x)} = {{\frac{c}{2\left( {c + \left( {x - b} \right)^{2}} \right)^{3/2}}{where}\quad b} = {{m - {1\quad{and}\quad c}} = {{3b} + {\frac{9}{4}.}}}}$The cumulative density function is given by${P_{g}\left( {x < y} \right)} = {\frac{1}{2}\left( {1 \pm \sqrt{1 - \frac{c}{c + \left( {y - b} \right)^{2}}}} \right)}$where the positive sign is taken if y>b and the negative sign otherwise.We can therefore draw randomly from g(x) by drawing u randomly from theinterval [0, 1], solving${u = {\frac{1}{2}\left( {1 \pm \sqrt{1 - \frac{c}{c + \left( {y - b} \right)^{2}}}} \right)\quad{for}\quad y}},$and making y be our random sample from g(x).

Now, let ${{h(x)} = \frac{g(x)}{f(x)}},$so that $\begin{matrix}{{h(x)} = \frac{\frac{c}{2\left( {c + \left( {x - b} \right)^{2}} \right)^{3/2}}}{\frac{x^{m - 1}e^{- x}}{\Gamma(m)}}} \\{= \frac{c\quad{\Gamma(m)}e^{x}}{2{x^{m - 1}\left( {c + \left( {x - b} \right)^{2}} \right)}^{3/2}}}\end{matrix}$and this is therefore minimised where$0 = {{h^{\prime}(x)} = {\frac{c\quad{\Gamma(m)}e^{x}}{2{x^{m - 1}\left( {c + \left( {x - b} \right)^{2}} \right)}^{3/2}}\left( {1 - \frac{m - 1}{x} - \frac{3\left( {x - b} \right)}{c + \left( {x - b} \right)^{2}}} \right)}}$which is where the final factor is zero, or$0 = {{{1 - \frac{m - 1}{x} - \frac{3\left( {x - b} \right)}{c + \left( {x - b} \right)^{2}}}\therefore{{x\left( {c + \left( {x - b} \right)^{2}} \right)} - {\left( {m - 1} \right)\left( {c + \left( {x - b} \right)^{2}} \right)} - {3\left( {x - b} \right)^{2}}}} = {{0\therefore{{\left( {x - b} \right)\left( {c + \left( {x - b} \right)^{2}} \right)} - {3\left( {x - b} \right)^{2}}}} = {{0\therefore{\left( {x - b} \right)\left( {c + \left( {x - b} \right)^{2} - {3\left( {x - b} \right)}} \right)}} = {{0\therefore{\left( {x - b} \right)\left( {\left( {x - b} \right)^{2} - {3\left( {x - b} \right)} + \left( {{3b} + \frac{9}{4}} \right)} \right)}} = 0}}}}$

The second factor has a negative discriminant since b>0, and cantherefore never be zero and is always positive; the only extremum ofh(x) is therefore at x=b and is a minimum; at that point we have${K = {{\min\limits_{x}\quad{h(x)}} = {{h(b)} = {\frac{c\quad{\Gamma(m)}e^{b}}{2{b^{m - 1}\left( {c + \left( {b - b} \right)^{2}} \right)}^{3/2}} = \frac{{\Gamma(m)}e^{b}}{2b^{b}c^{1/2}}}}}},$so the acceptance probability for the sample can be set at$\begin{matrix}{{p_{a}(x)} = {{{Kf}(x)}/{g(x)}}} \\{= {\frac{{\Gamma(m)}e^{b}}{2b^{b}c^{1/2}}\frac{2{x^{m - 1}\left( {c + \left( {x - b} \right)^{2}} \right)}^{3/2}}{c\quad{\Gamma(m)}e^{x}}}} \\{= {\frac{e^{b}}{b^{b}c^{3/2}}{x^{b}\left( {c + \left( {x - b} \right)^{2}} \right)}^{3/2}e^{- x}}}\end{matrix}$

Finally, we note that we will implement the acceptance decision bydrawing a uniformly distributed random number ν from the interval [0,1], and accepting our random sample if ν<P_(α)(x). The test for whetherν<p_(α)(x) or not can be speeded up in many instances by testing firstwhether ν is less than some other function of x that is simpler tocalculate and which is always less than p_(α)(x), and accepting thesample if it is.

m<1

In this case we set ${g(x)} = \left\{ \begin{matrix}\frac{{mx}^{m - 1}t}{t^{m}\left( {t + {me}^{- 1}} \right)} & {{{if}\quad 0} \leqq x \leqq t} \\\frac{{me}^{- x}}{t + {me}^{- 1}} & {{{if}\quad x} > t} \\0 & {otherwise}\end{matrix} \right.$where t=1−m.

The cumulative density function is${P_{g}\left( {x < y} \right)} = \left\{ \begin{matrix}\frac{y^{m}t}{t^{m}\left( {t + {me}^{- t}} \right)} & {{{if}\quad 0} \leq y \leq t} \\{1 - \frac{{me}^{- y}}{t + {me}^{- t}}} & {{{if}\quad y} > t} \\0 & {otherwise}\end{matrix} \right.$and we can therefore draw a random sample from g(x) by drawing auniformly distributed random number u from the interval [0, 1] andsolving P_(g)(x<y)=u for y.

Similarly to the derivation in the previous section, we can now set${h(x)} = {\frac{g(x)}{f(x)} = \left\{ \begin{matrix}\frac{{mt}\quad{\Gamma(m)}e^{x}}{t^{m}\left( {t + {me}^{- t}} \right)} & {{{if}\quad 0} \leq x \leq t} \\\frac{m\quad{\Gamma(m)}x^{t}}{t + {me}^{- t}} & {{{if}\quad x} > t}\end{matrix} \right.}$

The minima of both continuous sections are at the left hand end of thesection, i.e.${\underset{x \in I}{\arg\quad\min}\quad{g(x)}} = \left\{ \begin{matrix}0 & {{{if}\quad I} = \left\lbrack {0,t} \right\rbrack} \\{t +} & {{{if}\quad I} = \left( {x,t} \right\rbrack}\end{matrix} \right.$and are therefore at$K = \frac{{\Gamma(m)}{mt}}{t^{m}\left( {t + {me}^{- t}} \right)}$in both cases.

Therefore the acceptance probability${{p_{a}(x)}\quad{is}\quad e^{- y}\quad{if}\quad y} \leq {t\quad{or}\quad\left( \frac{y}{t} \right)^{m - 1}\quad{if}\quad y} > {t.}$

Finally, we note that we will implement the acceptance decision bydrawing a uniformly distributed random number ν from the interval [0,1], and accepting our random sample if ν<p_(α)(x). The test for whetherν<p_(α)(x) or not can be rearranged by first calculating$w = \left\{ \begin{matrix}x & {{{if}\quad x} \leq t} \\{t\quad{\log\left( \frac{x}{t} \right)}} & {{{if}\quad x} > t}\end{matrix} \right.$then testing whether ν<e^(−w). This test can then be speeded up in manyinstances by first testing whether $v > \frac{1}{1 + w}$(in which case we reject the sample and try all over again), or whetherν<1−w (in which case we accept the sample and have finished), and onlyif we reach this point do we then have to test whether ν<e^(−w).Random Samples from a Multinomial Distribution

We finally turn to efficient methods for sampling from a multinomialdistribution.

The multinomial distribution is analagous to the binomial distribution,but for tosses where there are more than 2 possible outcomes. Theparameters are N, the total number of tosses, and p, the vector ofprobabilities for the different possible results of a toss. Thedistribution is given by${P\left( {{n❘N},p} \right)} = {\frac{N!}{\prod\limits_{i = 1}^{I}\quad{n_{i}!}}{\prod\limits_{i = 1}^{I}\quad{p_{i}^{n_{i}}.}}}$Converting Production of a Single Multinomial Random Sample Into aSequence of Binomial Ones

The best way known to the author for sampling from this distribution isto convert the problem into a sequence of binomial random samples. Thefirst determines the values of n₁ and ${\sum\limits_{i = 2}^{I}n_{i}},$using as parameters for the binomial draw N and p₁. The kth determinesthe values of n_(k) and ${\sum\limits_{i = {k + 1}}^{I}n_{i}},$using as parameters for the binomial draw$N - {\sum\limits_{i = 1}^{k - 1}n_{i}}$and p_(k). That leaves the question of how to efficiently produce abinomial random sample.Production of Binomial Random Samples

The binomial distribution is the special case of the multinomialdistribution where I=2, except that usually (including here) the2-element parameter vector p is presented as the single element p=p₁(which uniquely determines p since the elements of the latter must sumto 1).

The best way known to the inventors for sampling from the binomialdistribution is to first order the values of$q_{n} = {\frac{N!}{\left( {N - n} \right)!}{p^{n}\left( {1 - p} \right)}^{N - n}}$in descending order (conceptually, without actually calculating themall), then to draw a uniformly distributed random number x from therange [0, 1], then to subtract from x in succession the values of q_(n)in succession starting with the largest. The output value n is the indexof the q_(n) which takes the remaining amount left in x below zero.

Implementation of this can take advantage of the fact that the methodworks almost as well if the q_(n) are used in a slightly different orderfrom the optimal one. We can therefore safely assume that${{\underset{n}{\arg\quad\max}\left( q_{n} \right)} = {{round}({Np})}},$and calculate this q_(n) first. Each subsequent one can then becalculated in order from either the one before or the one after it insuccession.

1. A method for determining one or more characteristics of a luminescentmaterial, the method comprising the steps of: a) irradiating theluminescent material with a pulse of excitation light, b) providing atriggering signal correlated to the pulse of excitation light; c)detecting with a photodetector a plurality of photons emitted from theluminescent material as result of the pulse of excitation light, thephotodetector providing an output signal upon detection of a photon; d)determining for each detected photon a photon arrival time and providingan output suitable for inputting to an analysing module wherein anoutput comprises zero, one, or more photon arrival times for eachexcitation; e) repeating steps a) to d) until a predetermined number ofexcitations has been performed; f) receiving said outputs in ananalysing module; and g) determining, in the analysing module,characteristic properties of the luminescent material by performing aprobabilistic analysis based on Bayesian inference.
 2. The method ofclaim 1, wherein the luminescent material has one or more luminescentmaterial characteristics such as luminescence lifetimes, luminescenceinitial intensities and ratios of luminescence initial intensities. 3.The method of claim 1, wherein the step of performing probabilisticanalysis based on Bayesian inference comprises determining from whichluminescent material each detected photon originates from.
 4. The methodof claim 1, wherein the step of performing statistical analysis based onBayesian inference comprises the steps of. i) setting a grid of valuesof luminescence lifetime α and initial intensity A; and ii) calculatingfor each value of α and A, a joint distribution, the joint distributionhaving the form${P\left( {A,\left. \alpha \middle| D \right.} \right)} = {A^{K + m_{A} - 1}\alpha^{m_{\alpha} - 1}{\mathbb{e}}^{- {A{({r_{A} + \frac{N}{\alpha}})}}}{\mathbb{e}}^{- {\alpha{({r_{\alpha} + {\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}t_{n,k}}}})}}}}$for determining of luminescence lifetime and initial intensity of asingle luminescent material.
 5. The method of claim 1, wherein the stepof performing probabilistic analysis based on Bayesian inferencecomprises including experimental artefacts in the Bayesian inference. 6.The method of claim 4, wherein the experimental artefacts comprises afirst timeperiod (δ) describing a PMT dead time and/or a secondtimeperiod (τ) describing blocking of a detector after an excitationand/or a third timeperiod (υ) describing the duration of recording aftereach excitation.
 7. The method of claim 1, wherein the step ofperforming probabilistic analysis based on Bayesian inference comprisestaking into account a first timeperiod (δ) describing a PMT dead time; asecond timeperiod (τ) describing blocking of a detector after anexcitation; and a third timeperiod (υ) describing the duration ofrecording after each excitation, and further comprises the steps of: i)setting a grid of values of luminescence lifetime α and initialintensity A; and ii) calculating for each value of α and A, a jointdistribution, the joint distribution having the form: $\begin{matrix}{{P\left( {A,\left. \alpha \middle| D \right.} \right)} \propto {A^{K + m_{A} - 1}\alpha^{m_{\alpha} - 1}}} \\{{\mathbb{e}}^{- {A{({r_{A} + {\frac{1}{\alpha}{({{N{({{\mathbb{e}}^{{- \alpha}\quad\tau} - {\mathbb{e}}^{{- \alpha}\quad\upsilon}})}} - {{({1 - {\mathbb{e}}^{{- \alpha}\quad\delta}})}{\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}{\mathbb{e}}^{{- \alpha}\quad t_{n,k}}}}}})}}})}}}} \\{{\mathbb{e}}^{- {\alpha{({r_{\alpha} + {\sum\limits_{n = 1}^{N}{\sum\limits_{k = 1}^{K_{n}}t_{n,k}}}})}}}}\end{matrix}$
 8. The method of claim 1, wherein steps a) to g) arerepeated until the statistical uncertainties of the measured propertiesare below predefined values or until a predetermined number ofrepetitions of steps a) to g) have been performed.
 9. The method ofclaim 1, wherein the step of performing probabilistic analysis based onBayesian inference comprises also inferring in which part of each timebin the true arrival time of each photon lies.
 10. A software programmodule for performing the method of claim
 1. 11. A measurement systemfor measuring luminescence characteristic properties of a luminescentmaterial, comprising a light source driving and controlling unit; and anexcitation light source for repeatedly irradiating the luminescentmaterial with pulses of excitation light, said light source driving andcontrolling unit providing a triggering signal correlated to each pulseof excitation light; at least one photodetector, for detecting aplurality of photons emitted from the luminescent material as result ofa pulse of excitation light, said photodetector providing an outputsignal upon the event of detection of a photon; means for photon timingdetermination coupled to receive said amplified output signal and saidtriggering signal, the means for photon timing determination determiningfor each emitted photon a photon arrival time and providing an output ofa plurality of photon arrival times, and means for analysing theplurality of photon arrival times, arranged to receive said plurality ofphoton arrival times outputted from said means from photon timingdetermination, and able to determine from the plurality of photonarrival times characteristic properties for the luminescent material,wherein the photon timing determination means determines and stores aplurality of arrival time measures corresponding to a plurality ofdetected photons, the plurality of photons resulting from the same pulseof excitation light; and wherein said means for analysing the pluralityof photon arrival times, utilises an algorithm for the determination ofcharacteristic luminescent properties of the luminescent material thatis based on Bayesian inference.
 12. The measurement system of claim 11,wherein said means for analysing the plurality of photon arrival timesis capable of performing probabilistic analysis based on Bayesianinference and inferring in which of each time bin the true arrival timeof each photon lies.
 13. The measurement system of claim 11, whereinsaid photon timing determination means is arranged to record the signalfrom the pre-amplifier at a sampling frequency above 1 GHz.